ad  685190 


AFFDL-TR-68-56 
VOLUME  I 


MAGIC:  AN  AUTOMATED  GENERAL  PURPOSE 
SYSTEM  FOR  STRUCTURAL  ANALYSIS 

VOLUME  I:  ENGINEER'S  MANUAL 

ROBERT  H.  MALLETT 
STEPHEN  JORDAN 

Bell  Aerosystems,  a  Textron  Company 


TECHNICAL  REPORT  AFFDL-TR-68-56 


JANUARY  1969 


n  D  C 
-3nn  nr?Irft 

APRS  «] 

G. 


This  document  has  been  approved  for  public 
release  and  sale;  its  distribution  is  unlimited. 


AIR  FORCE  FLIGHT  DYNAMICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO 


Reproduced  by  the 

CLEARINGHOUSE 
for  fcdor.il  Scientific  &  Technical 
Information  Springfield  Va  2215! 


This  Document  Contains 
Missing  Page/s  That  Are 
Unavailable  In  The 
Original  Document 


OR- 

ft  Icon  !C_- 

u\_i-  baLMC-' 


BEST 

AVAILABLE  COPY 


NOTICE 


When  Government  drawings,  specifications,  o-  other  data  are  used  for  any  purpose 
other  than  in  connection  with  a  definitely  relatec  Government  procurement  operation, 
the  United  States  Government  thereby  incurs  no  responsibility  nor  any  obligation 
whatsoever;  and  the  fact  that  the  Government  may  have  formulated,  furnished,  or  in 
any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not  to  be  regarded 
by  implication  or  otherwise  as  in  any  manner  licensing  the  holder  or  any  other  person 
or  corporation,  or  conveying  any  rights  or  permission  to  manufacture,  use,  or  sell  any 
patented  invention  that  may  in  any  way  be  related  thereto; 


This  document  has  been  approved  for  public  release  and  sale; 
its  distribution  is  unlimited. 


[SSiMfr 

'  J 

Icmn 

WHITE  SEC!lM  4  ' 

ilKT  SECTIM  □ 

ttMMWCn 

o 

jwnnwtiM  . 

1 . 

tat  . 

•«••••••«*••*  *  * 

j  aBTiiMTioa/iTAttMiun  coots  j 

ORT. 

1 

AVAIL  uA 

«r  «WAl 

Copies  of  this  report  should  not  be  returned  unless  return  is  required  by  security 
considerations,  contractual  obligations,  or  notice  on  a  specific  document. 


500  -  March  1969  -  C04S5  -  69-1521 


MAGIC:  AN  AUTOMATED  GENERAL  PURPOSE 
SYSTEM  FOR  STRUCTURAL  ANALYSIS 

VOLUME  I:  ENGINEER'S  MANUAL 

ROBERT  H.  MALLETT 
STEPHEN  JORDAN 


This  document  has  been  approved  for  public 
release  and  sale;  its  distribution  is  unlimited. 


FOREWORD 


This  report  was  prepared  by  Textron's  Bell  Aerosy stems,  Buffalo,  New  York, 
under  USAF  Contract  No.  AF  33  (615)-67-C-1505.  The  contract  was  initiated  under 
Project  No.  1467,  "Structural  Analysis  Methods,"  Task  No.  146702,  "Thermal  Elastic 
Analysis  Methods."  The  program  was  administered  by  the  Air  Force  Flight  Dynamics 
Laboratory  (AFFDL),  Air  Force  Systems  Command,  Wright- Patterson  Air  Force 
Base,  Ohio,  45433,  under  the  cognizance  of  Mr.  G.  E.  Maddux,  AFFDL  Program 
Manager.  The  program  was  carried  out  by  the  Structural  Systems  Department,  Bell 
Aerosystems,  during  the  period  15  March  1967  to  15  March  1968  under  the  direction 
of  Dr.  Robert  H.  Mallett,  Bell  Program  Manager. 

This  report,  "MAGIC:  An  Automated  General  Purpose  System  for  Structural 
Analysis,"  is  published  in  three  volumes;  "Volume  I:  Engineer's  Manual,"  Volume  II: 
User's  Manual,"  and  "Volume  III:  Programmer's  Manual."  The  manuscript  for 
Volume  I  was  released  by  the  authors  in  March  1968  ioi  publication. 

The  numerical  results  presented  in  this  report  were  obtained  at  the  Wright- 
Patterson  Air  Force  Base  Electronic  Data  Processing  Center.  The  utilization  of 
this  equipment  and  the  helpful  assistance  of  AFFDL  personnel  is  acknowledged. 

The  authors  wish  to  express  appreciation  to  colleagues  in  the  Advanced 
Structural  Design  Technology  Section  of  the  Structural  Systems  Department  for  their 
individually  significant,  and  collectively  indispensable,  contributions  to  this  effort. 
Special  acknowledgement  is  given  to  Mr.  Donald  Dupree  who  coordinated  the  final 
documentation  effort. 


The  authors  wish  to  express  appreciation  also  to  Miss  Beverly  J.  Dale  and 
Daniel  DeSantis,  for  the  expert  computer  programming  that  transformed  the  analyti¬ 
cal  development  into  a  practical  working  tool. 


This  technical  report  has  been  reviewed  and  is  approved. 


FRANCIS  J.  JANIK, 

Chief,  Theoretical  Mechanl 
Structures  Division 


Branch 


ABSTRACT 


An  automated  general  purpose  system  for  analysis  is  presented.  This  system, 
identified  by  the  acronym  "MAGIC"  for  "Matrix  Analysis  via  Generative  and  Inter¬ 
pretive  Computations,"  provides  a  flexible  framework  for  implementation  of  the 
finite  element  analysis  technology.  Powerful  capabilities  for  displacement,  stress 
and  stability  analyses  are  included  in  the  subject  MAGIC  System  for  structural 
analysis.  The  matrix  displacement  method  of  analysis  based  upon  finite  element 
idealization  is  employed  throughout.  Six  versatile  finite  elements  are  incorporated 
in  the  finite  element  library.  These  are:  frame,  shear  panel,  triangular  cross  sec¬ 
tion  ring,  toroidal  thin  shell  ring,  quadrilateral  thin  shell  and  triangular  thin  shell 
elements.  These  finite  element  representations  include  matrices  for  stiffness,  in¬ 
cremental  stiffness,  prestrain  load,  thermal  load,  distributed  mechanical  load  and 
stress.  The  MAGIC  System  for  structural  analysis  is  presented  as  an  integral  part 
of  the  overall  design  cycle.  Considerations  in  this  regard  include,  among  other 
things,  preprinted  input  data  forms,  automated  data  generation,  data  confirmation 
features,  restart  options,  automated  output  data  reduction  and  readable  output  dis¬ 
plays.  Documentation  of  the  MAGIC  System  is  presented  in  three  parts:  namely. 
Volume  I:  Engineer's  Manual;  Volume  II:  User's  Manual;  and  Volume  HI:  Pro¬ 
grammer's  Manual.  Volume  I  is  the  primary  technical  document.  Included  are  a 
general  technical  discussion  of  the  MAGIC  System,  an  outline  of  the  theoretical 
framework,  statement  of  the  individual  finite  element  representations,  and  illustra¬ 
tive  analyses  for  evaluation  of  each  finite  element  representation.  Volume  n  contains 
instructions  for  the  preparation  of  input  data  and  interpretation  of  output  data  with 
examples  drawn  from  the  illustrations  presented  in  Volume  I.  Volume  IE  is  designed 
to  facilitate  implementation,  operation,  modification  and  extension  of  the  MAGIC 
System. 
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1.  INTRODUCTION 


Bell  Aerosystems  has  been  active  in  the  development  of  automated  structural 
analysis  tools  based  upon  the  finite  element  technology  since  the  late  1950's.  In  a  con¬ 
tractual  outgrowth  of  this  internal  development  activity,  Bell  furnished  a  series  of  com¬ 
puter  programs  to  the  Air  Force  Flight  Dynamics  Laboratory  (AFFDL)  in  1963.  These 
programs,  described  in  References  1  through  6,  became  an  integral  part  of  structural 
analysis  practices  at  AFFDL  and  at  numerous  other  recipient  governmental  and  pri¬ 
vate  organizations. 

Advances  in  computer  software  and  hardware  signaled  the  impending  obselescence 
of  the  foregoing  computer  programs  for  structural  analysis  in  1966.  Attempts  to  sal¬ 
vage  these  programs  by  direct  modifications  to  the  coding  proved  discouraging.  More¬ 
over,  newly  established  technological  advances  strongly  recommended  development  of 
a  second  generation  finite  element  capability  for  structural  analysis. 

In  the  light  of  the  situation  just  described,  Bell  undertook,  in  March  of  1967,  to 
implement  an  advanced  general  purpose  system  for  Matrix  Analysis  via  Generative 
and  Interpretive  Computations  (MAGIC)  at  AFFDL.  This  MAGIC  System  for  struc¬ 
tural  analysis,  described  herein,  was  planned  to  provide,  as  a  minimum,  the  capability 
of  the  prior  set  of  Bell  computer  programs.  The  capability  ultimately  built  into  the 
MAGIC  System  is  actually  far  more  powerful  than  the  former  programs  taken  collec¬ 
tively.  For  example,  structures  characterized  by  "on  the  order  of"  2000  degrees-of- 
freedomcanbe  accommodated  in  contrast  to  tire  former  500  degrees-of-freedom  limit. 

Documentation  of  the  MAGIC  System  for  structural  analysis  is  presented  in 
three  volumes.  The  subject  volume  (Volume  1)  is  the  primary  technical  report.  The 
major  sections  of  this  report  are  described  in  the  following  paragraphs.  Separate 
supplementary  volumes  are  provided  to  facilitate  utilization  of  the  MAGIC  System. 
Volume  n,  the  User’s  Manual^,  includes  detailed  specifications  for  the  preparation 
of  input  data,  along  with  illustrative  examples.  Volume  in,  the  Programmer's  Manual^®), 
presents  information  on  the  organization  of  the  computer  program  as  well  as  its 
operational  characteristics. 

A  general  description  of  the  MAGIC  System  for  structural  analysis  is  included 
in  Section  2.  Particular  attention  is  given  to  definition  of  the  overall  organization  of 
the  system.  A  key  element  of  this  organization  is  seen  to  be  the  versatile,  AFFDL 
sponsored,  FORTRAN  Matrix  Abstraction  Technique  (FORMAT  II)  described  in  Refer¬ 
ences  9  through  12.  Emphasis  is  also  given  in  this  section  to  special  data  management 
features  which  facilitate  efficient  utilization  of  the  MAGIC  System  such  as  preprinted 
input  data  forms. 

Section  3  of  this  primary  technical  report  outlines  the  theoretical  bases  employed 
in  derivation  of  the  finite  element  representations  and  in  development  of  the  analysis 
procedures.  A  total  of  six  finite  elements  ore  incorporated  in  the  Element  Library  of 
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the  MAGIC  System;  namely,  frame,  shear  panel,  triangular  cross  section  ring,  toroi¬ 
dal  thin  shell  ring,  quadrilateral  thin  shell  and  triangular  thin  shell  elements.  The 
computational  procedures  outlined  in  Section  3  include  displacement,  stress  and 
stability  analyses. 

Sections  4  through  9  present  statements  of  the  matrices  which  comprise  the  in¬ 
dividual  finite  element  representations.  In  general;  stiffness,  incremental  stiffness, 
pressure  load,  thermal  load,  and  stress  matrices  are  provided.  Sections  4  through  9 
also  include  numerical  evaluations  of  the  respective  finite  elements.  These  evalua¬ 
tions  take  the  form  of  series  of  selected  example  problems. 

The  body  of  this  technical  report  is  concluded  with  a  general  retrospective  dis¬ 
cussion  in  Section  10..  The  MAGIC  System  is  given  critical  review.  Limitations  are 
discussed  and  guidelines  for  utilization  are  presented. 
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2,  TECHNICAL  DISCUSSION 


A.  INTRODUCTION 

Automated  general  purpose  capabilities  promise  to  revolutionize  analysis  and 
design  practices.  The  matrix  methods  of  analysis  based  upon  discrete  element  ideal¬ 
ization  provide  the  suitable  theoretical  basis.  High  speed  data  processing  devices 
establish  the  economic  feasibility.  Powerful  automated  tools  for  analysis  and  design 
have  already  been  derived  from  these  resources.  Experience  accumulated  in  the  de¬ 
velopment  and  application  of  these  tools  has  evolved  a  conceptual  framework  suitable 
for  generalization. 

Expansions  which  traverse  traditional  boundaries  between  the  specialized  disci¬ 
plines  of  mechanics,  improvements  which  provide  firm  theoretical  bases  for  consist¬ 
ent  mathematical  models^13),  and  extensions  which  automate  design  iterations^14)  are 
now  well  defined.  New  data  management  concepts  which  facilitate  data  handling^-®), 
matrix  abstraction  instructions  which  simplify  programming^®),  and  hardware  devices 
which  enable  convenient  display  (and  communication)^1®)  have  also  emerged.  The  Im¬ 
plementation  of  all  these  generalizations  within  the  framework  of  realistic  hardware 
design  poses  a  stimulating  challenge. 

The  advanced  general  purpose  system  for  Matrix  Analysis  via  Generative  and 
Interpretive  Computations  (MAGIC)  which  is  described  herein  was  developed  in  ac¬ 
ceptance  of  the  foregoing  challenge.  This  MAGIC  System  furnishes  the  specific 
structural  analyses  capability  sought  and,  at  the  same  time,  provides  a  versatile  con¬ 
ceptual  framework  to  facilitate  the  foregoing  generalizations.  Accordingly,  general 
concepts  are  given  consideration  in  the  following  discussion  as  well  as  specific  fea¬ 
tures  of  the  MAGIC  System  for  structural  analysis. 

B.  ANALYSIS  TECHNOLOGY 


l'ne  finite  element  approach  to  structural  analysis  is  consistently  stated  in 


Section  3  within  the  framework  of  the  variational  methods  of  continuum  mechanics. 
Within  this  framework,  discretization  can  be  referenced  to  zones  designed  to  facilitate 
the  construction  of  admissible  displacement  function  mode  shapes.  Elementary  illus¬ 
trative  physical  models  arising  from  such  an  idealization  into  zones  are  shown  in 
Figures  la  and  lc.  Admissible  assumed  displacement  functions  written  individu¬ 
ally  for  each  zone,  when  taken  collectively,  form  admissible  assumed  displacement 
functions  whose  field  of  definition  is  the  entire  structure. 


These  physical  models,  formed  by  subdivision  into  zones,  may  be  equivalently 
viewed  as  assemblies  of  discrete  structural  elements  interconnected  such  that  ap¬ 
propriate  interclement  continuity  is  maintained.  For  example,  the  physical  models 
shown  in  Figures  la  and  lc  may  be  equivalently  viewed  as  assemblies  of  the  discrete 
structural  elements  of  Figures  lb  and  Id,  respectively.  It  is  this  latter  viewpoint, 
taken  herein,  which  makes  evident  the  generality  of  the  finite  element  methods  of 
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Figure  1.  Illustrative  Structural  Idealizations 


analysis.  Idealization  into  zones  or  structural  components  enables  the  systematic 
treatment  of  large  scale  complex  structures  as  assemblages  of  large  numbers  of 
common  elementary  structural  components. 

Mathematical  models  are  formulated  for  selected  elementary  structural  com¬ 
ponents  of  parametrically  specified  shape.  These  finite  element  representations  are 
then  given  specific  dimensions  to  form  building  blocks  appropriate  to  structures  of  a 
general  problem  class.  Interconnection  of  adjacent  elements  is  provided  for  by  the 
construction  of  displacement  function  mode  shapes  with  gridpoint  displacement  func¬ 
tion  quantities  as  undetermined  coefficients.  Taking  these  gridpoint  displacement 
degrees-of-freedom  common  to  adjacent  elements  establishes  their  interconnection. 

The  referencing  of  the  structural  idealization  to  elementary  physical  components 
leads  naturally  to  specification  of  descriptive  data  with  respect  to  these  individual 
elements.  Variations  in  dimensions  such  as  thickness  are  accommodated  by  the 
specification  of  distinct  values  for  individual  elements.  Material  property  variations 
arising  from  lamination  or  temperature  degradation  are  accommodated  by  element 
related  characterizations  of  materials. 

Distributed  loadings  are  also  processed  by  element  in  order  to  account  for 
variations.  Elementary  distributions  are  assumed  over  individual  elements  in  much 
the  same  way  that  displacement  function  mode  shapes  are  constructed.  Intensities  of 
distributed  loadings  such  as  pressure,  temperatures  and  prestrain  are  prescribed  at 
gridpoints.  These  intensities  are  then  transformed  into  work  equivalent  forces  via 
the  assumed  distributions. 

The  foregoing  comments  have  indicated  the  facility  with  which  finite  element 
idealization  accommodates  problematical  variations  in  geometry,  material  and  applied 
loading.  It  is  useful  to  emphasize  this  point  further  by  examination  of  the  overall 
computational  process. 

The  basic  computational  flow  of  a  finite  element  stress  analysis  is  illustrated 
in  Figure  2.  The  important  feature  to  be  noted  in  this  flow  chart  is  that  the  mathe¬ 
matical  description  of  a  structural  system  (Block  2)  is  generated  independently  of  the 
construction  of  the  objective  mathematical  model  for  the  structural  system  (Block  3). 
That  is,  physical  description  (elastic  constants,  pressures,  etc.)  is  referenced  to  the 
individual  zones  or  finite  elements  and  transformed  to  appropriate  element  mathemat¬ 
ical  representation  without  regard  to  total  structure  configuration  and  boundary  con¬ 
ditions.  It  is  primarily  this  separation  which  accounts  for  the  generality  of  the  dis¬ 
crete  element  method  in  regard  to  both  complexity  and  broad  applicability. 

Regarding  complexity,  referencing  of  problem  description  to  individual  discrete 
elements  enables  convenient  consideration  of  variations  in  geometry,  sizing  dimen¬ 
sions,  material  properties,  applied  loadings,  and  boundary  conditions.  Regarding  ap¬ 
plicability,  this  is  limited  only  by  the  suitability  of  the  discrete  elements  made  avail¬ 
able  for  idealization. 
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Figure  2.  Structural  Analysis  Computational  Flow 


A  variational  point  of  view  is  maintained  throughout  the  subject  analysis  proc¬ 
ess.  Specifically,  the  principle  of  potential  energy  is  employed.  The  well  known 
Rayleigh- Ritz  assumed  mode  method  of  analysis  is  invoked  to  generate  the  desired 
algebraic  expressions  for  the  element  energy  functions.  Then,  these  are  summed  to 
obtain  the  energy  function  for  the  total  structure.  The  objective  governing  equations 
follow  immediately  by  executing  the  variation  of  the  total  energy.  The  principal  ad¬ 
vantage  of  maintaining  the  variational  viewpoint  throughout  this  process  is  that  the 
matrices  involved  enjoy  explicit  and  complete  labeling  at  every  step.  The  theoretical 
framework  is  outlined  in  detail  in  Section  3.  Therein,  the  analysis  processes  are 
given  explicit  definition  in  terms  of  matrices. 

C.  FINITE  ELEMENTS 

The  MAGIC  System  incorporates  the  six  finite  elements  shown  in  Figure  3; 
namely,  frame,  shear  panel,  triangular  cross  section  ring,  toroidal  thin  shell  ring, 
quadrilateral  thin  shell  and  triangular  thin  shell  elements.  These  elements,  taken 
collectively,  enable  the  idealization  of  most  structures. 

The  set  of  matrices  embodied  in  each  element  representation  determines  the 
type  of  analyses  which  can  be  performed.  In  the  MAGIC  System,  a  complete  element 
representation  is  taken  to  include  matrices  for  stiffness,  incremental  stiffness,  pres¬ 
sure  load,  prestrain  load,  thermal  load  and  stress.  Moreover,  provision  has  been 
made  for  additional  element  matrices  such  as  consistent  mass  matrices. 

The  frame  element  is  a  conventional  "beam  theory"  finite  element.  This  ele¬ 
ment  is  well  suited  to  the  idealization  of  planar  and  space  frames.  An  eccentric  con¬ 
nection  feature  is  incorporated  in  this  frame  element  representation  to  facilitate 
utilization  as  a  shell  stiffener  element.  The  frame  element  is  also  appropriate  to 
planar  and  space  trusses. 

The  truss  specialization  of  the  frame  element  is  particularly  useful  in  combi¬ 
nation  with  the  quadrilateral  shear  panel  element.  The  quadrilateral  shear  panel  ele¬ 
ment  simulates  the  action  of  a  thin  panel  in  diagonal  tension.  The  effective  extension- 
al  stiffness  is  allocated  to  truss  elements.  Such  axial  force  member-shear  panel 
idealizations  have  found  extensive  application  in  the  analysis  of  airframe  structures. 

The  triangular  cross  section  ring  element  is  one  of  *he  earliest  and  best  known 
finite  element  models.  This  versatile  element  enables  realistic  idealization  of  thick- 
walled  axisymmetric  structures  of  arbitrary  profile. 

The  representation  of  the  triangular  cross  section  ring  incorporated  in  the 
MAGIC  System  is  basically  the  same  as  the  original  mortel^7)  although  several  useful 
generalizations  have  been  introduced.  One  of  these  is  the  orthotropic  material  cap¬ 
ability  with  data  specified  orientation  of  material  axes. 


The  integrations  conducted  in  formulating  this  element  also  serve  to  set  it  apart. 
A  precise  integration  is  carried  out  when  the  radial  dimension  of  the  cross  section  is 
not  small  relative  to  the  diameter  of  the  ring.  vVhen  the  radial  dimension  of  the  ring 
is  small  relative  to  the  diameter,  an  approximate  integration  is  carried  out  in  accord 
with  that  in  the  conventional  representation. 

The  MAGIC  System  representation  of  the  triangular  cross  section  ring  embodies 
matrices  for  pressure  and  prestrain  load  as  well  as  for  stiffness  and  stress.  A  par¬ 
ticularization  of  the  prestrain  load  vector  is  included  to  facilitate  thermal  loading. 

The  thin  shell  elements  incorporated  in  the  MAGIC  System  are  particularly 
noteworthy  since  they  have  not  been  presented  previously  in  the  open  technical  litera¬ 
ture.  The  toroidal  thin  shell  ring  represents  a  substantial  improvement  over  the 
predecessor  conic  thin  shell  ring^®).  In  contrast  to  the  latter,  the  toroidal  ring  yields 
accurate  predictions  of  stresses  for  relatively  coarse  idealizations.  In  applications 
where  the  double  curvature  of  the  toroidal  ring  is  not  required,  it  specializes  to  conic 
and  cylindrical  configurations.  Moreover,  the  toroidal  ring  reduces  easily  to  a  cap  or 
end  closure  element. 

The  quadrilateral  and  triangular  sets  of  thin  shell  elements  incorporated  in  the 
MAGIC  System  provide  an  unprecedented  capability  for  the  analysis  of  thin  membrane, 
plate  and  shell  structures.  The  arbitrary  shape  of  these  elements  enables  efficient 
idealization  of  complex  configurations  and  gridwork  refinement.  Supplementary  mid¬ 
side  gridpoints  are  optionally  available  to  facilitate  local  gridwork  refinement. 

Interelement  continuity  is  assured  between  elements  of  common  and  companion 
type.  As  a  consequence,  recourse  to  convergence  criteria  is  often  permitted.  The 
variation  in  strains  built  into  these  elements  yields  accurate  stress  predictions  rela¬ 
tive  to  predecessor  elements^. 

Many  additional  special  features  are  included  in  this  set  of  thin  shell  elements 
and  in  the  other  elements  as  well;  for  example,  arbitrary  material  axes,  arbitrary 
stress  axes,  plane  strain  option,  etc.  It  is  features  such  as  these  which  establish  the 
MAGIC  System  as  a  practical  analysis  tool  as  opposed  to  simply  a  large  scale  finite 
element  computer  program. 

A  separate  section  of  this  report  is  devoted  to  the  presentation  of  each  of  these 
finite  element  representations.  The  reader  is  directed  to  the  introductions  within 
these  sections  for  further  description  of  the  finite  elements  and  their  representations. 

D.  PROGRAMMING  TECHNOLOGY 

Useful  insight  into  the  nature  of  the  finite  element  based  MAGIC  System  for 
stress  analysis  can  be  gained  from  examination  of  the  conceptual  organizational  chart 
shown  in  Figure  4.  This  chart  illustrates  the  modularization  which  is  fundamental  to 
general  purpose  program  organization.  Overall  efficiency  is  achieved  by  this  modular¬ 
ization  in  much  the  same  way  that  complex  electronic,  mechanical,  and  even  structural 
systems  are  modularized  to  maximize  versatility  and  maintainability. 
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Figure  4.  General  Purpose  Program  Organizational  Chart 


The  Libraries  shown  in  Figure  4  are  particularly  noteworthy.  These  represent 
a  higher  level  in  a  hierarchy  of  modularization  in  that  they  build  in  entire  series  of 
optional  modules.  The  simultaneous  availability  of  alternatives,  achieved  by  standard¬ 
ized  module  interfaces,  prov'des  numerous  benefits.  The  most  obvious  benefits  are 
those  derived  from  flexibility.  Standardization  also  gives  a  repetitiveness  to  program 
developmental  phases  that  enhances  efficiency  and  reliability.  Engineering  interfaces 
reflect  this  standardization  to  advantage  as  well.  These  factors  contribute  importantly 
to  the  favorable  cost  effectiveness  of  the  MAGIC  System  which  is  discussed  in  Section 
11. 


The  conceptual  organization  shown  in  Figure  4  reaches  slightly  beyond  that  of 
the  subject  MAGIC  System,  This  enables  a  more  comprehensive  discussion  of  the 
relevant  programming  technology  and  gives  perspective  to  the  actual  organization  of 
the  MAGIC  System.  Variances  between  the  organization  of  Figure  4  and  that  of  the 
MAGIC  System  are  delineated  clearly. 

The  nature  and  function  of  the  individual  program  modules  are  described  in  the 
following  paragraphs.  These  descriptions  also  indicate  the  position  of  the  modules  in 
the  logical  flow  of  an  analysis.  As  a  consequence,  the  individual  module  descriptions, 
taken  collectively,  yield  the  objective  overall  picture  of  the  programming  technology 
intrinsic  to  the  MAGIC  System. 

1.  Resident  Operating  System 

The  Resident  Operating  System  controls  and  coordinates  job  processing.  It 
normally  contains  such  subsystems  as  input/output  routines,  external  storage  super¬ 
visors,  language  compilers  and  assemblers  and  system  accounting  routines.  Example 
Resident  Operating  Systems  are:  IBSYS  for  IBM  7090  and  7094,  OS  for  the  IBM 
System/360,  EXEC  for  the  UNIVAC  1108  and  SCOPE  for  the  CDC  6600. 

Machine  compatibility  has  been  insured  by  the  exclusive  use  of  FORTRAN 
IV  in  the  MAGIC  System.  The  absence  of  machine  or  assembler  language  from  every 
portion  of  the  program  eliminates  most  problems  of  machine  dependency  and  imple¬ 
mentation  difficulties.  Thus,  even  though  the  program  is  a  system  in  itself,  it  is  de¬ 
signed  to  function  under  the  control  of  the  normal  operating  system  resident  on  a 
machine. 


Avoidance  of  machine  dependency  also  prevents  optimum  utilization  of  aux¬ 
iliary  direct  access  storage  units.  However,  the  overall  organization  of  Figure  4  is 
designed  to  accommodate  this  generalization.  The  conceptual  logic  implied  in  the 
chart  implr  ?  the  addition  of  more  modules  under  the  Executive  Monitor  than  just  the 
FORMAT  Monitor.  In  this  way,  direct  access  dependency  could  be  incorporated  into  a 
monitor  on  the  same  level  as  the  FORMAT  Monitor.  Other  matrix  control  systems 
could  also  be  placed  under  the  Executive  Monitor.  Every  addition  would  also  auto¬ 
matically  inherit  the  capabilities  of  the  underlying  modules. 
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2.  Executive  Monitor 


The  Executive  Monitor  is  the  highest  level  of  control  within  the  MAGIC 
System.  This  module  controls  location  of  the  Problem,  Execution,  and  Material 
Libraries.  In  addition,  the  Executive  Monitor  has  sole  control  over  maintaining  and 
accessing  the  Execution  Library.  The  primary  function  of  the  Executive  Monitor  is  to 
coordinate  the  Libraries  in  conjunction  with  selection  of  the  appropriate  submonitor 
as  directed  by  the  application.  Since,  at  the  present  time,  the  FORMAT  Monitor  is  the 
only  submonitor,  the  alternative  Libraries  are  placed  under  its  control  and  the  Execu¬ 
tion  Monitor  is  not  required. 

3.  Problem  Library 

The  Problem  Library  takes  the  form  of  a  magnetic  tape  prepared  for  the 
Analyst.  Multiple  problems  are  accumulated  on  the  Problem  Library  tape.  An  entry 
in  the  Problem  Library  includes  a  complete  record  of  the  input  data  specification, 
selected  intermediate  results,  and  the  output  data  specification.  Control  of  the  tape  as 
regards  access  for  subsequent  additions,  deletions,  calculation,  or  displays  resides 
with  the  Analyst.  The  Problem  Library  serves  as  a  flexible  interface  between  pro¬ 
gram  and  Analyst  thereby  providing  an  opportunity  for  effective  data  management. 

Input  data  sets  are  made  self  generating  insofar  as  is  possible.  This  is 
effected  in  a  preprocessing  phase  and  the  complete  data  set  is  recorded  in  tl  3  Problem 
Library.  After  approval  of  the  input,  the  Analyst  can  invoke  the  Problem  Library  to 
continue  the  execution.  Placement  of  intermediate  results  in  the  Problem  Library 
provides,  in  the  same  way,  for  economical  recovery  at  certain  milestones  in  the 
solution  process. 

Localized  design  modifications  and  gridwork  refinements  can  be  accom¬ 
modated  without  dealing  with  the  entire  structure.  This  is  particularly  important 
where  multiple  thermal  loads  are  considered.  Partitioning  can  sometimes  be  designed 
to  circumvent  ill  conditioning.  Very  large,  very  sparse  matrices  are  avoided,  as  are 
long  continuous  executions.  Generally  it  can  be  said  that  the  Problem  Library,  partic¬ 
ularly  when  coupled  with  substruoturing,  allows  analysis  operations  to  be  broken  into 
manageable  units. 

In  summary,  the  Problem  Library  (Block  3)  of  Figure  4,  is  accessible  from 
any  module  below  it  to  obtain  data  from  previous  problems  and  store  data  for  the  cur¬ 
rent  problem.  A  Problem  Library  may  be  generated  in  the  present  MAGIC  System  to 
the  extent  provided  by  the  availability  of  restart  points  in  the  analysis  process  as  de¬ 
scribed  under  the  FORMAT  Monitor.  No  provision  is  built  in  to  conduct  analyses  by 
substructures. 
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4.  Execution  Library 


The  Execution  Library  is  designed  to  build  in  alternative  abstraction  in¬ 
struction  sequences.  Entries  in  this  module  represent  procedures  such  as  displace¬ 
ment,  stress  and  stability  analyses.  In  addition  to  standard  built-in  analyses,  non¬ 
standard  analyses  can  be  conducted  simply  by  defining  an  entry  in  the  Execution 
Library.  Revision  or  deletion  of  this  module  is  controlled  by  the  Executive  Monitor. 

The  broad  variety  of  analyses  encountered  in  practice  actually  embody 
relatively  few  computations  which  are  unique;  rather,  an  extensive  commonality 
exists.  It  is  this  commonality  that  enables  the  efficient  development  and  operation  of 
automated  capabilities  which  are  general  purpose  in  the  sense  of  multiple  types  of 
analyses. 


This  module  is  not  built  into  the  present  MAGIC  System.  Abstraction  in¬ 
struction  sequences  are  included  in  the  input  data  deck  to  effect  the  desired  type  of 
analysis. 


5.  Material  Library 

The  specification  of  mechanical  and  physical  material  properties  can  be  a 
burdensome  task.  This  is  particularly  true  in  the  case  of  laminated  materials  or  in 
the  presence  of  thermal  degradation  of  material  properties.  Accordingly,  the  Material 
Library  is  a  very  useful  feature  of  the  MAGIC  System.  This  Material  Library  is  sim¬ 
ple  in  concept;  yet,  its  availability  can  save  time  measurable  in  man-days  against  a 
single  problem.  In  contending  with  design  changes  and  multiple  thermal  load  condi¬ 
tions,  the  Materia]  Library  is  virtually  indispensible. 

The  Material  Library  takes  the  form  of  a  magnetic  tape  which  is  a  perma¬ 
nent  data  set  available  for  interrogation  by  the  MAGIC  System.  The  Executive  Monitor 
is  the  natural  control  level  for  additions,  modifications  and  deletions  to  the  Material 
Library.  In  the  absence  of  the  Executive  Monitor,  this  function  is  served  by  the 
Structural  System  Monitor  in  the  present  MAGIC  System.  Updating  of  the  Material 
Library  may  be  conducted  as  a  separate  execution  or  as  an  integral  part  of  the 
analysis  process. 

A  complete  set  of  temperature  referenced  properties  for  a  material  con¬ 
stitutes  an  entry  in  the  Material  Library.  Each  entry  in  the  Material  Library  is  taken 
to  include  material  designation,  lock  code,  elastic  constants,  coefficients  of  thermal 
expansion  and  mass  density.  Provision  is  made  for  data  at  up  to  nine  temperature 
levels.  Linear  interpolation  is  employed  in  interrogation  of  the  Material  Library  for 
material  property  values  at  a  specified  temperature  level.  Material  anisotropy  is 
assumed  as  well  as  temperature  dependence. 
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6.  FORMAT  Monitor 


In  the  absence  of  an  Executive  Monitor,  the  control  functions  and  respons¬ 
ibilities  of  the  Executive  Monitor  are  handled  by  or  delegated  to  the  FORMAT  Monitor. 
In  addition,  the  FORMAT  Monitor  carries  out  its  normal  functions. 

The  FORMAT  Monitor  controls  the  selection  and  usage  of  the  underlying 
modules  within  the  confines  permitted  by  the  Execution  Monitor.  At  each  transfer 
point  between  the  underlying  modules  the  FORMAT  Monitor  will  make  a  logical  de¬ 
cision,  based  upon  information  returned  from  the  module,  regarding  the  continuance 
or  discontinuance  of  processing.  Termination  of  processing  is  determined  voluntarily 
by  the  Analyst  unless  unrecoverable  error  conditions  are  encountered  by  a  module. 

The  FORMAT  Monitor  contains  the  correlation  table  between  external 
storage  devices  and  their  respective  FORMAT  functions.  The  Analyst  has  at  his 
command  the  option  to  revise  the  correlation  table  for  any  given  application.  The 
FORMAT  Monitor  has  the  assignment  of  processing  any  such  revisions. 

Restart  capabilities  are  also  controlled  by  the  FORMAT  Monitor  as  directed 
by  the  Analyst.  By  generating  the  desired  abstraction  instruction  sequence  and  re¬ 
questing  pertinent  Information  to  be  saved,  the  Analyst  has  at  his  command  flexible 
restart  capabilities.  For  example,  in  the  contexts  of  a  structral  system,  element 
matrices  may  be  generated,  saved  and  the  problem  restarted  at  a  later  level  of  anal¬ 
ysis.  Another  example  of  restart  would  be  to  utilize  the  option  of  termination  after 
the  Structural  System  Input  Data  has  been  read  and  interpreted.  Saving  of  this  inter¬ 
preted  input  would  allow  the  Analyst  to  examine  the  input  printout  and  restart  the 
problem  without  the  necessity  of  reinserting  the  original  data  for  reading  and  inter¬ 
pretation. 


Operating  under  the  FORMAT  Monitor,  the  basic  computational  flow  of  the 
program  starts  at  the  Preprocessor  Monitor,  passes  the  Execution  Monitor  and  then 
to  the  Structural  System  Monitor  which  ends  the  cycle  by  returning  control  to  the 
Execution  Monitor.  In  this  way,  multiple  data  decks  may  be  batched  in  a  single 
MAGIC  System  execution. 

7.  Preprocessor  Monitor 

The  Preprocessor  Monitor  interprets  problem  specification  data  pertinent 
to  program  setup.  The  processing  involved  includes  specification  of  (a)  master  input 
tapes,  (b)  master  output  tapes,  (c)  analysis  header  labels,  (d)  problem  header  labels 
and  (e)  page  size  for  printout.  Matrices  provided  via  input  data  cards  are  read  and 
stored  within  the  Preprocessor  Monitor.  Other  functions  of  the  Preprocessor  Monitor 
are  accomplished  through  its  three  underlying  modules. 


8.  Abstraction  Instruction  Compiler 

The  Abstraction  Instruction  Compiler  interprets  the  abstraction  instructions 
and  extracts  matrix  names,  operation  codes,  scalars  and  statement  numbers  in  the 
process.  These  quantities  are  stored  in  packed  form  and  returned  to  the  Preprocessor 
Monitor  for  use  by  the  Instruction  Logic  Supervisor.  Serious  compilation  errors  may 
terminate  execution  at  this  point. 

9.  Machine  Resources  Allocator 

The  Machine  Resources  Allocator  partitions  the  available  internal  storage 
into  a  program  area  and  work  area.  This  module  also  assigns  program  functions  to 
the  external  storage  facilities  available.  The  four  possible  program  functions  for  an 
external  storage  device  are  instruction  storage,  master  input  unit,  master  output  unit 
and  input/output  utility  unit.  These  allocations  of  storage  areas  are  based  upon  pro¬ 
gram  and  application  requirements.  If  no  master  input  or  master  output  units  are 
needed,  their  function  reverts  to  input/output  utility. 

10.  Instruction  Logic  Supervisor 

The  Instruction  Logic  Supervisor  scans  the  information  assembled  by  the 
Abstraction  Instruction  Compiler  and  the  Machine  Resources  Allocator  and  creates  a 
logical  path  for  the  Execution  Monitor.  At  the  same  time  an  optimum  external  storage 
assignment  is  made  for  each  matrix  named  in  an  abstraction  instruction  in  the  know¬ 
ledge  of  the  logical  path  to  be  followed.  The  Instruction  Logic  Supervisor  takes  into 
account,  in  this  process,  such  consideration  as  the  channel  addresses  of  external 
storage  facilities,  number  of  external  storage  facilities,  capacities  of  external  storage 
facilities  and  combinations  of  input  and  output  matrices  of  abstraction  instructions. 

The  result  is  an  optimum  utilization  of  available  machine  resources  for  the  sequence 
of  operations  released  to  the  Execution  Monitor  from  the  Preprocessor  Monitor. 

11.  Execution  Monitor 

The  Execution  Monitor  follows  the  path  specified  by  the  Preprocessor  Moni¬ 
tor  accessing  the  underlying  modules  to  perform  the  prescribed  operations.  Operations 
can  be  performed  on  matrices  up  to  the  order  2000.  The  efficient  utilization  of  machine 
storage  resources  is  assured  by  the  setup  passed  from  the  Preprocessor  Monitor. 

The  Execution  Monitor  will  terminate  processing  if  any  of  the  rules  of 
matrix  algebra  are  violated.  Matrices  are  stored  by  columns  complete  with  matrix 
name,  dimensions  and  sign.  If  a  column  of  a  matrix  is  less  than  50%  dense,  it  is 
stored  in  compressed  format.  The  modular  form  allows  ease  of  insertion  of  additional 
matrix  manipulative  or  generative  operations. 
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12.  Algebraic  Matrix  Operation 

The  Algebraic  Matrix  Operation  module  is  essentially  a  library  of  routines 
for  matrix  manipulation.  This  library  includes  routines  for  addition,  subtraction, 
multiplication,  transposed  multiplication,  scalar  multiplication,  transposition,  inver¬ 
sion,  equation  solving  by  elimination,  equation  solving  by  iteration,  and  eigenvalue/ 
eigenvector  extraction. 

Each  of  the  above  operations  is  incorporated  into  a  separate  module  and 
all  except  the  eigenvalue  operation  have  out  of  core  capability. 

13.  Nonstandard  Matrix  Operation 

The  Nonstandard  Matrix  Operation  is  essentially  a  library  of  routines  to 
effect  nonstandard  matrix  manipulation.  Included  in  this  library  are  routines  to  raise 
each  element  within  a  matrix  to  a  specified  power,  locate  maximum  or  minimum 
values  in  a  row  or  column,  adjoin  two  matrices  column  wise,  and  multiply  two  matrices 
element  by  element. 

14.  Special  Function  Modules 

The  Special  Function  Modules  constitute  a  library  of  routines  to  effect  non- 
algebraic  operations.  Included  in  this  library  are  routines  to  print,  skip  ahead  upon 
encountering  a  null  matrix,  select  the  best  condition  columns  from  a  triangular 
matrix,  and  solve  the  selected  set  of  simultaneous  equations,  form  a  diagonal  matrix 
from  a  row  or  column  matrix,  and  rename  a  matrix. 

15.  Structural  System  Monitor 

The  Structural  System  Monitor  is  the  matrix  generator  of  the  structural 
analysis  capability  provided  by  the  MAGIC  System.  Machine  storage  resources  are 
allocated  to  this  module  by  the  Preprocessor  Monitor.  Matrices  describing  a  struc¬ 
tural  system  are  released  from  this  monitor  for  the  conduct  of  the  matrix  manipula¬ 
tion  phase  of  the  structural  analysis  process.  The  Structural  System  Monitor  together 
with  its  underlying  modules  comprise  the  major  portion  of  the  MAGIC  System  for 
structural  analysis. 

16.  Structure  Data  Preprocessor 

The  Structure  Data  Preprocessor  is  the  principal  input  data  interface  be¬ 
tween  the  MAGIC  System  and  the  Structural  Analyst.  As  such,  the  nature  of  this 
module  is  described  in  the  subsection  E,  ’'Program/Analyst  Interfaces." 

The  basic  function  of  the  Structure  Data  Preprocessor  is  to  read  and  inter¬ 
pret  all  data  describing  the  idealized  structral  model  and  to  make  this  data  available 
for  the  generation  of  structural  matrices  via  the  Structural  System  Monitor.  The 
interpretation  function  carried  out  by  the  Structure  Data  Preprocessor  is  substantial 
since  data  sets  are  designed  to  be  internally  generated  insofar  as  is  possible. 
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An  optional  execution  interruption  is  provided  at  completion  of  the  structural 
data  preprocessing.  The  completed  set  of  structural  data  is  printed  for  examination  by 
the  Analyst.  Then,  upon  approval  of  the  input  data,  the  analysis  process  is  restarted. 

17.  Utility  Library 

The  Utility  Library  is  an  elementary  interpretive  system  in  the  form  of  a 
collection  of  FORTRAN  subroutines.  Computational  routines  which  are  common  to 
several  element  matrix  generation  procedures  are  placed  in  the  Utility  Library  to 
avoid  a  duplication  of  programming.  An  extensive  commonality  exists  among  the 
generation  procedures  even  for  diverse  types  of  discrete  elements.  Exploitation  of 
this  commonality  via  the  Utility  Library  contributes  measurably  to  the  efficient  de¬ 
velopment  of  the  Element  Library  in  the  MAGIC  System.  Included  in  the  Utility 
Library  are  routines  for  numerical  integration,  interpolation,  specialized  structural 
print  and  algebraic  operations  for  small  size  matrices. 

18.  Element  Library 

The  Element  Library  is  the  heart  of  the  MAGIC  System  for  structural  anal¬ 
ysis.  Each  entry  in  this  library  represents  a  finite  element  model.  A  call  on  the  Ele¬ 
ment  Library  causes  numerical  generation  of  certain  matrices  of  a  complete  element 
representation. 

The  availability  in  the  Element  Library  of  suitable  elements  for  idealization 
determines  the  applicability  of  an  analysis  system  to  different  classes  of  structure. 
Moreover,  the  set  of  matrices  embodied  in  each  element  representation  determines 
the  type  of  analyses  which  can  be  performed.  In  the  absence  of  versatile  Element 
Libraries,  even  the  best  matrix  and  tape  interpretive  systems  yield  sterile  analysis 
capabilities. 

The  six  finite  element  models  incorporated  in  the  Element  Library  of  the 
MAGIC  System  and  the  set  of  element  matrices  provided  were  described  in  the  pre¬ 
ceding  subsection  C.  Experience  has  shown  this  Element  Library  to  provide  a  power¬ 
ful  capability  for  structural  analysis. 

E.  PROGRAM/ANALYST  INTERFACES 

Discussion  of  the  MAGIC  System  L.r  structural  analysis  is  not  complete  without 
some  comment  on  the  program/Analyst  interfaces.  The  acceptance  of  automated 
analysis  tools  by  stress  analysts  hinges  importantly  on  the  simplicity  of  these  inter¬ 
faces.  The  first  interface  encountered  by  the  Analyst  is  with  the  Preprocessor  Moni¬ 
tor.  The  basic  instruction  sequence  to  be  executed  passes  through  this  interface  from 
the  Analyst  to  the  program.  These  instructions  consist  of  a  sequence  of  mathematical 
equations  to  be  performed.  An  abstraction  instruction  sequence  for  linear  stress 
analysis  is  illustrated  in  Figure  5.  Such  instruction  sequences  may  be  constructed  at 
the  volition  of  the  Analyst  and  executed  to  perform  a  wide  variety  of  computations. 
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DISPLACEMENT  AND  STRESS  ANALYSIS  INSTRUCTION  SEQUENCE 

MATL8A, LOADS, TR , TA, KEL , FEL , SEL , SZALEL ,  ,  ,  *  =  .U5EF.Q4. 

PRINT  OLTPUT  MATRICES 

PRINT  (O.O.F.'COND.  ,E6,I  LOADS 
PRINT  ( REDCOF ,C .0 . F . , E6, )  TR 
PRINT  iNSYS  ,N0kSUM,E6f )  TA 
PRINT  I  ROM  »C0L  ,6,,»  KEL 

PRINT  (ROM  'COL  ,E6,I  FEL 

PRINT  (ROM  'COL  ,E6,)  SEL 

PRINT  (ROM  tCOL  ,Eb,>  SZALEL 

FORM  TAR  MATRIX  ( ASSEMBLY  AND  APPLICATION  OF  BOUNDARY  COND^ » 

TRT  *  TR  .VRANSP. 

TAR  *  TA  .TMULT.  TRT 

ASSEMBLE  ANC  REDUCE  ELEMENT  STIFFNESS  MATRICES 

KTEMP  *  KEL  . TMULT.  TAR 
STIFF  *  TAR  .TMULT.  KTEMP 
PRINT  (FORCE  «CISP.  ,  «)  STIFF 

ASSEMBLE  ANC  REDUCE  ELEMENT  APPLIED  LOADS 

FTELAR  *  TAR  .TMULT.  FEL 
PRINT  ( REDDOF ,COND.  ,  ,)  FTELAR 

APPLY  BOUNDARY  CONDITIONS  TO  SYSTEM  LOADS 

LOADR  *  TR  .MULT.  LOAOS 
PRINT  (REOOQF'CONO.  ,  *)  LOADR 

COMBINE  ELEMENT  AND  SYSTEM  LOAOS 

TLOAO  »  FTELAR  .jADO.  LOADR 
PRINT  (REOCOF.CONO."  TLOAO 

SOLVE  FCR  DISPLACEMENTS 

DISPR  »  STIFF  .SEQEL.  TLOAD 
PRINT  ( REDCOF, CONO.  ,  ?  I  DISPR 

SOLVE  FCR  ELEMENT  STRESSES 

STREL  *  SEL  .MLLT.  TA* 

STRESF  »  STREL  .MULT.  DISPR 

STRESS  «  STRESF  .SUBT.  SZALEL 
PRINT  (NRSEL  ,COND.  ,  , )  STRESS 

SOLVE  FOR  ELEMENT  FORCES 

FQRCEL  «  KTEMP  .MULT.  DISPR 
FORCES  «  FORCEL  .SUBT.  FEL 
PRINT  (O.O.F.'COND.  »  , )  FORCES 

SOLVE  FGR  SYSTEM  REACTIONS 

REACTN  »  TA  .MILT.  FORCES 
REACT  a  REACTN  .SUBV.  LOADS 
PRINT  (O.O.F.'COND.  ,  ,)  REACT 

Figure  5.  Illustrative  Abstraction  Instruction  Listing 
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BELLOOIO 

BELL  0020 

BELL  0030 

BELL  0040 

BELL0050 

BELL0060 

BELL0070 

BELL0080 

BELL 0090 

BELL 0100 ‘ 

BELL0110 

BELL0120 

BELL0130 

BELL0140 

BELL  01 50 

BELL  0160 

BELL0170 

BELLOl'80 

BELL0190 

8ELL0200 

BELL  0210 

’BELL  0220' 
BELL0230 
BELL0240 
BELL0250 
BELL0260 
BELL0270 

’ BELL 0280" 
BELL0290 
BELL 0360 
BELL0310 
BELL 0320 
BELL0330 
BELL0345  " 
BELL0350 
BEU.0360 
BELL0370 
BELL 03 80 
BELL0390 

'bellTWFo 

BELL0410 
BELL0420 
BELL0430 
BELL 0440 
8ELL0450 
BELL 04 60 
BELL0470 
BELL0430 
8ELL0490 
BELL0500 

BELL0510 
BELL0520 
BELL0530 
BELL9540 
BELL0550 
BELL0560" 
BELL0570 
Bell  0586 
BELL0590 
SELLO&OO 
BELL061C 
BELLO  620 
8EL10630 
8ELL0640 


Executions  may  be  terminated  and  restarted  at  the  corresponding  exit  and  entry 
points  of  any  abstraction  instruction.  Input  data,  intermediate  results  or  final  results 
can  be  automatically  saved  in  this  way.  Then,  with  the  retrieval  of  this  data,  comput- 
ation  can  be  resumed. 

The  second  program/Analyst  interface  encountered  is  with  the  Structural 
System  Monitor.  This  is  the  primary  input  data  interface  of  the  MAGIC  System  for 
structural  analysis.  Experience  has  shown  that  significant  portions  of  the  labor  and 
computer  costs  of  analyses  are  occasioned  by  incomplete  or  improper  specification  of 
problem  input  data.  In  recognition  of  this,  special  features  are  associated  with  the 
MAGIC  System  to  facilitate  the  confirmation  of  problem  data  prior  to  execution.  In¬ 
cluded  are  annotated  input  forms,  data  consistency  checks,  and  an  option  to  read,  com¬ 
plete  and  write  the  input  data  prior  to  attempting  execution. 

Preprinted  input  data  forms  are  essential  to  the  reliable  specification  of  data. 
These  forms  provide  a  labeled  entry  position  for  all  data  items  which  gives  engineer¬ 
ing  definition  to  the  quantities  requested.  Control  options  are  selected  simply  by  a 
mark  (X).  These  provisions  help  to  minimize  occurrences  of  incomplete  specifications 
of  problem  data. 

The  printed  input  forms  take  advantage  of  a  special  MODAL  data  card  feature. 

The  MODAL  card  feature  enables  data -prescribed  initialization  of  tables.  Explicit 
data  requirements  are  thereby  limited  to  specification  of  exceptions  to  the  MODAL 
initialization. 

In  addition  to  the  MODAL  card,  a  data  Repeat  option  is  available.  When  utilized, 
data  from  the  previous  point  is  retained  for  the  indicated  point.  The  combination  of 
the  MODAL  card  and  the  Repeat  option  significantly  reduces  the  volume  and  complexity 
of  input. 

The  input  forms  also  embody  permanent  label  cards  which  automatically  precede 
subsets  of  data,  thereby  allowing  flexibility  in  the  arrangement  of  the  subsets  of  data 
to  form  the  total  input  data  deck.  Data  associated  with  options  not  exercised  are 
simply  omitted.  This  is  particularly  useful  when  a  problem  is  bemg  restarted  at  an 
advanced  stage  of  computation. 

A  data  confirmation  preprocessing  phase,  with  problem  execution  suppressed, 
is  a  recommended  practice  in  utilization  of  the  MAGIC  System.  In  this  data  process¬ 
ing  execution,  explicit  data  is  read  and  Implied  data  is  generated.  For  example, 
MODAL  card  completions  are  conducted  and  material  properties  are  interpolated 
from  the  Material  Library.  Consistency  of  all  the  data  is  checked  ar.d  a  complete 
record  of  the  data  is  recorded  for  restart  and  printed  for  inspection. 

There  are  basically  two  types  of  output  provided  by  the  MAGIC  System.  The 
first  is  matrix  print  provided  from  the  Special  Function  Print  Module.  This  encom¬ 
passes  all  output  external  to  the  Structural  System  Monitor.  A  standard  format  is 
employed  to  print  matrices. 
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The  second  basic  f.ype  of  output  is  that  provided  from  within  the  Structural 
System  Monitor.  Output  from  this  module  includes  a  list  of  the  completed  input  data 
with  self-explanatory  engineering  labels.  In  addition,  intermediate  results  employed 
in  checkout  are  optionally  available  in  the  completed  program. 

F.  SIZE  CHARACTERISTICS 

The  size  characteristics  of  the  MAGIC  System  are  twofold:  first,  there  are  the 
size  characteristics  of  the  program  itself  and  second,  those  associated  with  the  prob¬ 
lem  solving  capability.  Considering  the  former,  the  MAGIC  System  contains  212  sub¬ 
routines  (approximately  25,000  FORTRAN  IV  source  cards)  logically  designed  into  89 
overlay  links  on  an  IBM  7090  with  32,000  words  of  storage.  The  overlay  design  re¬ 
flects  the  optimum  use  of  available  storage  yet  maintains  respectable  execution  effi¬ 
ciency. 

The  MAGIC  System  offers  large  scale  capability  with  no  penalties  to  small 
applications  due  to  the  fact  that  out  of  core  operations  are  not  utilized  unless  the  mag¬ 
nitude  of  the  application  requires  them.  The  size  of  the  program  has  necessitated  use 
of  SUBSYS,  a  package  which  improves  the  loading  capabilities  of  IBSYS,  on  the  7090/ 
94.  In  addition  to  allowing  the  program  to  be  loaded,  SUBSYS  allows  the  program 
overlay  load  tape  to  be  saved,  thereby  improving  execution  time.  Also,  SUBSYS  allows 
programs  to  be  executed  back  to  back  without  passing  through  the  IBLDR  section  of 
IBJOB  for  each  program.  On  the  7090  under  SUBSYS  the  program  is  actually  divided 
into  three  segments:  Preprocessor.  Execution  and  Structural  System.  Third  genera¬ 
tion  computers,  such  as  System/360  and  UNIVAC  1108,  have  the  capabilities  of 
SUBSYS  incorporated  into  their  resident  operating  system. 

The  scale  of  the  analysis  capability  provided  via  the  MAGIC  System  can  be 
characterized  as  "on  the  order  of"  2000  displacement  degrees-of-freedom.  Other 
relevant  maximum  size  characteristics  are  1000  discrete  elements,  1000  grid  points 
and  10  applied  load  conditions.  Matrices  which  r.re  card  input  may  be  of  order  2000 
x  2000  and  contain  up  to  4500  single  precision  real  non-zero  elements  on  a  32,000  word 
machine. 

The  MAGIC  System  needs  a  minimum  of  eight  external  storage  units  to  operate, 
distributed  into  the  following  functions:  one  unit  assigned  as  Instruction  storage  for 
the  Execution  Monitor,  one  unit  assigned  as  a  Master  Input  Unit,  one  unit  assigned  as 
a  Master  Output  Unit,  and  five  units  assigned  as  Input/Output  Utility  Units.  Every 
effort  should  be  made  to  make  the  most  external  storage  units  possible  available, 
since  any  increase  in  the  available  storage  units  increases  execution  efficiency. 

The  stated  maximum  size  characteristics  apply  to  the  linear  stress  analysis 
capability  of  the  MAGIC  System.  A  stability  analysis  capability  is  also  included  in  the 
MAGIC  System,  as  with  the  linear  stress  analysis,  and  explicit  matrix  statement  of 
the  stability  analysis  procedure  is  given  in  Section  3.  The  number  of  displacement 
degrees-of-freedom  which  can  be  accommodated  in  the  eigenvalue  stability  analysis  is 
limited  to  130.  The  other  size  characteristics  stated  for  the  stress  analysis  remain 
applicable. 


G.  SPECIAL  FEATURES 


Many  features  have  been  built  into  the  MAGIC  System  which  are  not  fundamental 
to  a  finite  element  computer  program  but  which  are  essential  to  a  general  purpose 
analysis  system  for  practical  structures.  Foremost  among  these  features  is  a  great 
variety  of  transformation  matrices.  Material  axes  transformations  are  provided  to 
accommodate  arbitrary  axes  of  orthotropy.  Stress  axes  transformations  enable  the 
referencing  of  output  displays  to  convenient  axis  systems.  Grid  point  axes  transform¬ 
ations  account  for  irregular  boundary  conditions  and  allow  pseudo-curvilinear  dis¬ 
placement  variables.  Eccentric  connection  transformations  provide  for  realistic 
modeling  of  frame  joints  and  shell  stiffeners.  Finally,  grid  point  suppression  trans¬ 
formations  are  included  to  eliminate  unwanted  element  grid  points  prior  to  assembly. 

A  second  feature  of  special  interest  is  the  element  repeat  feature.  There  are 
actually  two  levels  of  element  repeat.  The  first  is  a  repeat  of  element  data.  Under 
this  option,  all  calculations  proceed  as  usual,  but  the  repeated  provision  c ?  identical 
element  extra  data  cards  is  avoided.  The  second  level  of  element  repeat  is  element 
matrix  repeat  and  this  is  the  more  powerful  option  by  far.  Under  this  option,  the  ele¬ 
ment  matrices  of  the  prior  element  are  simply  carried  forward  as  those  of  the  present 
element;  no  calculation  is  carried  out.  ulearly  under  this  option,  a  great  saving  in 
input  data  specification  is  realized  and  important  savings  in  calculation  can  be  realized 
as  well.  The  extent  to  which  the  input  data  can  be  reduced  by  the  element  matrix  re¬ 
peat  feature  is  made  clear  in  the  User's  Manual. 

A  useful  element  load  condition  scalar  is  associated  with  the  multiple  load  con¬ 
dition  capability  of  the  MAGIC  System.  Element  load  conditions  arise  in  load  condition 
number  one.  A  multiplicative  constant  is  then  data  prescribed  for  all  subsequent  load 
conditions.  This  scalar  controls  the  participation  of  the  element  loading.  With  this 
feature,  a  total  load  system  can  be  decomposed  into  several  parts  and  behavior  pre¬ 
dictions  can  be  obtained  conveniently  against  these  as  well  as  against  the  total  load 
system.  This  feature  is  particularly  useful  in  separating  effects  of  thermal  and 
mechanical  applied  load  combinations. 

The  majority  of  the  special  features  embodied  in  the  MAGIC  System  are  explained 
best  within  a  specific  context.  Accordingly,  with  the  exception  of  the  few  included  here 
for  special  emphasis,  such  features  are  treated  as  an  integral  part  of  other  report 
sections.  Many  are  disclosed  in  Volume  II  of  this  report  in  the  process  of  explaining 
items  of  input  data  and  interpreting  example  problem  output  data. 
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3.  THEORETICAL  FRAMEWORK 


A.  INTRODUCTION 

The  matrix  methods;  of  analysis  based  upon  discrete  element  idealization  have 
been  the  subject  of  an  extensive  body  of  technical  literature  and,  more  recently, 
entire  books  as  well.  (19>  20).  This  documentation  obviates  the  need  for  detailed 
theoretical  development  herein.  Nevertheless,  in  the  interest  of  clarity  and  com¬ 
pleteness,  presentation  of  the  discrete  element  representations  incorpoi'ated  within 
the  MAGIC  System  is  prefaced  in  this  Section  by  general  symbolic  statement  of  the 
analysis  processes.  This  gives  explicit  definition  to  the  methodology  and  notation 
employed. 

Statement  of  the  analysis  processes  is  separated  into  three  parts,  Firstly, 
consideration  is  given  to  the  discrete  element  representations.  Then,  having  given 
definition  to  the  discrete  element  matrices  employed,  the  steps  executed  by  the 
MAGIC  System  in  the  conduct  of  a  linear  stress  analysis  are  described.  Lastly,  the 
stability  analysis  process,  which  is  an  extension  of  the  linear  stress  analysis,  is 
presented. 

B.  DISCRETE  ELEMENT  MATRICES 

1.  Fundamental  Requirements 

The  development  of  a  discrete  element  representation  is  essentially  a 
problem  in  elasticity.  Accordingly,  the  fundamental  requirements  to  be  satisfied 
are  those  of: 

(a)  Equilibrium, 

(b)  Material  Behavior, 

(c)  Compatibility,  and 

(d)  Boundary  Conditions. 

It  is  convenient  to  approach  the  satisfaction  of  these  requirements  for  a  discrete 
element  variationally  by  way  of  the  principle  of  potential  energy^*)  which  states 
that: 


Of  all  possible  displacement  states  within 
as  given  admissible  class  {  8  }  » that  which 
makes  the  total  potential  energy  <J>  p  j  g  j 
stationary,  satisfies  the  equilibrium  re-J 
quirements  and  is  t  ^  actual  displacement 
state  {s}  *,  i.e. 
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Furthermore,  if 


a$P  {&} 

"Isj 


=  0  j  -  1,2,  ...,  n. 

{8}" {8}* 


%  ({»}*)<*,({»}) 


(1) 


(2) 


for  all  {8}  in  some  neighborhood  of  ,  then  associated  equilibrium  position 
is  stable. 

2.  Discretization 

The  foregoing  statement  of  the  principle  of  potential  energy  is  expressed 
in  terms  of  a  finite  number  of  displacement  variables  /  g  1  implying  prior  discreti¬ 
zation  of  the  potential  energy  functional.  Discretization  of  the  potential  energy 
functional  is  effected  in  accordance  with  the  well  known  Rayleigh--Ritz  techniques 
by  the  introduction  of  assumed  displacement  mode  shapes.  Admissibility  conditions 
must  be  imposed  on  the  characteristics  of  these  displacement  mode  shapes  to  assure 
satisfaction  of  certain  fundamental  requirements. 

The  fundamental  requirement  of  compatibility  of  strains  is  provided  for 
subsequently  in  this  development  by  expression  of  the  strains  in  terms  of  displace¬ 
ments.  Since  the  functional  dependence  of  strains  upon  displacements  involves 
differentiation,  continuity  requirements  arise  as  criteria  of  admissibility  to  be 
satisfied  in  the  construction  of  displacement  mode  shapes.  It  should  be  emphasized 
that  these  continuity  requirements  remain  applicable  across  discrete  element 
boundaries  ^ 22 K 

The  foregoing  interelement  continuity  admissibility  conditions  are  peculiar 
to  the  discrete  element  method  of  analysis.  The  admissibility  requirements  associated 
with  conventional  applications  of  the  Rayleigh-Ritz  techniques  apply  as  well.  The 
definition  of  general  systematic  procedures  for  constructing  displacement  functions 
within  the  collective  confines  of  these  fundamental  requirement  related  admissibility 
conditions  has  proved  to  be  an  elusive  goal.  However,  significant  progress  in  this 
direction  has  been  made  by  the  use  of  unconventional  and  curvel inear  coordinate 
systems  and  interpolation  formuiae(23,  24,  25). 

Practical  considerations  involved  in  the  selection  of  assumed  displacement 
functions  go  beyond  the  problem  of  admissibility.  Of  particular  importance  is  the 
number  of  displacement  degree3-of- freedom  to  be  associated  with  an  element.  The 
provision  of  degree- of- freedom  in  excess  of  the  number  required  to  establish  ad¬ 
missibility  is  attractive  in  that  it  reduces  the  number  of  elements  required  in  ideali¬ 
zations  in  order  to  maintain  a  certain  level  of  precision  and  correspondingly  reduces 
the  input  data  preparation. 
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Improvement  in  stress  predictions  is  also  realized  as  a  consequence  of 
including  additional  degrees-of-freedom  in  an  element  representation.  Furthermore, 
it  has  beer,  demonstrated  on  certain  example  problems  that  improved  predictions  of 
displacement  behavior  can  be  obtained  with  fewer  total  degrees-of-freedom  if  the 
number  of  degrees-of-freedom  associated  with  an  individual  element  are  increased^26)- 

These  attractive  advantages  of  higher  order  assumed  displacement  functions 
are  achieved  at  the  expense  of  simplicity,  which  has  been  a  primary  recommendation 
of  the  discrete  element  methods.  This  characteristic  has  been  somewhat  obscured  by 
the  trend  toward  advanced  geometrically  complex  discrete  elements  pursued  in  the 
interest  of  eliminating  structure  idealization  errors.  The  additional  increment  in 
complexity  of  mode  definition,  formulation,  checkout,  specification,  and  numerical 
expression  introduced  by  extra  element  degrees-of-freedom  severely  handicaps 
attempts  to  achieve  the  aforementioned  advantages. 

As  a  final  comment  regarding  criteria  for  selection  of  assumed  displacement 
functions  it  is  pertinent  to  note  that  many  practical  structures  have  obvious  physical 
definition  in  terms  of  panels  and  stiffners.  A  lesser  element  gridwork  would  require 
prohibitively  complex,  problem  orientated,  stiffened  panel  discrete  elements.  At  the 
same  time  the  increase  in  accuracy  afforded  by  a  higher  order  panel  element  repre¬ 
sentation  is  unwarranted  in  most  problems  of  this  type.  Thus,  it  is  concluded  that 
the  most  significant  advancements  in  element  representations  will  continue  to  stem 
from  elimination  of  structure  idealization  error  rather  than  reduction  of  element  dis¬ 
cretization  error. 

The  actual  process  of  constructing  displacement  mode  shapes  begins  with 
the  definition  of  a  convenient  set  of  coordinate  axes  for  the  discrete  element  model. 
Then,  the  boundaries  of  the  element  are  given  parametric  description.  Polynomial 
mode  shapes  are  the  type  customarily  chosen  to  represent  the  displacement  functions 
within  the  parametrically  described  boundaries  of  a  discrete  element.  With  reference 
to  the  selected  element  coordinate  axes,  such  assumed  displacement  functions  can  be 
written  symbolically  as 


{“<>}  -  [0<>H4 


(3) 


where 


u  j  is  the  vector  of  displacement  functions, 

B  is  the  matrix  of  mode  shapes,  and 
#  •  is  the  vector  of  mode  shape  participation  coefficients. 


The  participation  coefficients  |^3jin  the  assumed  displacement  modes  are  referred 
to  as  "field  coordinate"  displacement  degrees-of-freedom.  These  field  coordinates 
are  commonly  retained  throughout  the  algebraic  development  of  a  discrete  element 
representation;  however,  in  order  to  effect  assembly  of  elements  (establish  interelement 


continuity)  it  is  necessary  to  transform  to  gridpoint  displacement  degrees-of-freedom 
/ 8  j,  .  This  transformation  results  from  a  straightforward  application  of  inter¬ 
polation  theory.  The  displacement  functions  are  particularized  to  the  selected  grid- 
point  quantities  j  8  j  thereby  yielding, 


{ 8s  J =  [  r8/s]  t  ^ 


(4) 


The  objective  transformation  is  then  obtained  by  the  inversion  of  this  relation,  i.e. 


{*}•  [i>8]  {sg} 


(5) 


The  gridpoint  displacement  degrees -of-freedom{s}  are  generally  defined 
with  respect  to  coordinate  systems  on  the  individual  discrete  elements.  Frequently, 
a  number  of  further  displacement  coordinate  transformations  are  then  necessary  to 
obtain  degrees-of-freedom  which  are  suitable  for  assembly  and  convenient  for  inter¬ 
pretation.  All  such  transformations  are  given  explicit  definition  within  the  individual 
discrete  element  representations;  however,  two  are  common  to  most  elements  and 
are  described  here. 


Generally,  it  is  necessary  to  transform  to  a  global  Cartesian  set  of  co¬ 
ordinate  axes.  This  system,  common  to  all  discrete  elements  of  an  idealized 
structure,  is  suitable  for  interconnection  of  the  elements.  The  transformation  re¬ 
lation  to  obtain  gridpoint  displacement  degrees-of-freedom  1 8  S1  referenced  to  global 
axes  takes  the  form 


{Sg}  [  rgs  ]{  8s  } 


(6) 


in  which  the  transformation  matrix  L^gsJ  consists  of  submatrices  of  direction 
cosines. 


Boundary  conditions  on  displacement  quantities  not  aligned  with  the  global 
axes  require  special  point-related  coordinate  axes  for  these  gridpoints.  Taking  the 
associated  coordinate  axes  transformation  for  a  gridpoint  as, 

{‘.HJ.toj  <7) 

the  transformation  to  gridpoint  axis  displacement  degrees-of-freedom  is  given  by 

{ss  }=[rsq  ]{sq} 

Transformations  of  this  type  are  employed  simply  to  facilitate  interpretation  of 
the  results  in  many  cases. 
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It  is  useful  to  conclude  comment  on  the  construction  of  displacement 
function  mode  shapes  by  collecting  the  foregoing  transformations.  The  result  is 


where 


(10) 


Customarily,  the  formulative  process  ;'.s  carried  forward  using  the 
field  coordinate  displacement  degrees-of-freedom  ■{  }  and  then  Equation  9  is  invoked 

to  obtain  the  discrete  element  matrices  with  respect  to  the  gridpoint  displacement 
degrees-of-freedom  {S  q}*  The  matrices  which  actually  participate  in  this  collective 
transformation  ^  ]  vary  from  element  to  element. 

3.  Equilibrium 


The  pi'inciple  of  potential  energy  was  introduced  to  facilitate  satisfaction 
of  the  fundamental  requirements  for  a  discrete  element.  Having  examined  the 
nature  of  the  discretization  implied  in  the  statement  of  the  energy  principle, 
attention  is  returned  to  assuring  satisfaction  of  these  fundamental  requirements. 

It  is  clear  from  the  statement  of  the  principle  of  potential  energy  that  this 
variational  approach  circumvents  explicit  consideration  of  equilibrium  requirements. 
The  equilibrium  requirements  arise  naturally  in  the  Euler  equations  of  the  variation 
process.  This  is  an  important  advantage  of  the  method. 

4.  Material  Behavior 


Preceeding  in  the  order  listed  at  the  outset,  the  second  fundamental  re¬ 
quirement  to  be  satisfied  in  the  elasticity  problem  posed  by  a  discrete  element  is 
that  of  material  behavior.  Linear  elastic  behavior,  governed  by  a  generalized  Hooke's 
law,  is  assumed,  i.e. 


{{<<>}  -{*!<>}} 


where 


is  the  stress  state, 
is  the  state  of  strain 

is  the  elastic  property  characterization,  and 
is  the  prestrain  state. 


(11) 
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In  recognition  of  the  increasing  utilization  of  high  performance  particulate 
and  fibrous  composite  materials,  material  anisotropy  is  provided  for  in  defining  this 
stress-strain  relation.  The  availability  of  material  property  data  generally  limits 
material  specifications  to  orthotropic,  at  most.  However,  the  application  of  a  rotational 
transformation  in  order  to  reference  the  material  characterization  to  the  geometric 
coordinate  axes  of  a  discrete  element  tends  to  fill  the  material  property  matrices. 

For  this  reason,  no  terms  in  elastic  [  E]  and  thermal  {  a}  property  characterization 
matrices  are  assumed  zero. 

5 .  Compatibility 

Satisfaction  of  the  third  fundamental  requirement,  compatibility,  is  provided 
for  by  expressing  strains  in  terms  of  the  displacements.  Interpretation  of  this  require¬ 
ment  in  terms  of  admissibility  conditions  on  displacement  mode  shapes  was  discussed 
previously  and  appropriate  functions  are  assumed  available  at  this  point.  The  intro¬ 
duction  of  these  displacement  mode  shapes  (Equation  3)  into  the  relevant  strain-dis¬ 
placement  equations  enables  expression  of  the  strains  in  terms  of  the  discrete  element 
field  coordinate  displacement  degrees- of- freedom,  i.e. 

{e  < )}  =  [c  (>]{#}  <“> 


Nonlinear  terms  have  been  omitted  in  this  set  of  strain-displacement 
relations.  These  will  be  given  special  consideration  subsequently. 

6 .  Boundary  Conditions 

The  final  fundamental  requirements  which  must  be  established  are  the 
boundary  conditions.  Force  boundaries  need  not  be  given  explicit  consideration 
since  these  are  accommodated  implicitly  by  the  variational  process.  Displacement 
boundary  conditions,  on  the  other  hand,  must  be  imposed.  Expression  of  the  element 
displacement  mode  shapes  in  terms  of  boundary  displacement  provides  for  the  simple 
imposition  of  these  boundary  conditions. 


7.  Potential  Energy 


Proceeding  toward  algebraic  expression  of  a  discrete  element  representa¬ 
tion,  it  remains  to  give  definition  to  the  potential  energy  function.  The  strain  energy 
density,  dU,  which  is  basic  to  the  potential  energy,  is  defined  as 


l  d'J{<r} 


(13) 


Invoking  the  relation  governing  material  behavior  obtain  expression  of  the 
total  strain  energy  as, 
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(14) 


v-f  (tI«][eH«H*  J[eH  })dv 

*  \r 


Substituting  the  relation  governing  compatibility  (Equation  12)  obtain,  in  the  pre¬ 
sence  of  distributed  mechanical  loading,  |  P  ]•  the  total  potential  energy  function  in 

*p  '/v  "7  l^J[c()f[E][c()]{^} 

~[p\  [B(  }]T{p)  {15) 


This  substitution  of  the  assumed  displacement  functions  into  the  element 
total  potential  energy  functional  and  the  subsequent  integration  over  the  volume 
comprise  a  major  part  of  the  effort  associated  with  the  derivation  of  a  discrete 
element  representation.  The  procedure  is  conceptually  simple  though  algebraically 
complex.  Indicating  the  integration  symbolically,  obtain  an  algebriac  expression  for 
the  element  total  potential  energy  as 


-l  *>]{%}  ™ 

where  /*  ~  „  t 

[k]-/v  [  Ct)]  [  E][  Ct)]dV 

(17) 

{~&l  [c<>]T[E]{ei  (,}av 

(18) 

{%}  =  /jB<  )} 

(19) 

These  element  m.  trices  in  the  potential  energy  expression  are  referenced 
to  the  field  coordinate  displacement  degrees-of- freedom.  The  previously  defined 
transformation  (Equation  9)  is  introduced  to  obtain  the  element  matrices  in  the 
potential  energy  expression  with  reference  to  selected  gridpoint  displacement 
degrees- of- freedom,  i.e. 

V  i  l8JIKH#,  H8,  JfoH8, 

J{  FP  }  <2C> 
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where 


{  F«H  r/3q  ]  T{?«  } 
{'pMrPq]T{?P} 


At  this  point  the  objective  matrices  governing  behavior  of  a  discrete 
element  follow  immediately  by  executing  the  variation  of  the  potential  energy 
function,  i.e. 

MK  Hr«HFp}‘  lFc} 

where 

is  the  element  stiffness  matrix, 
is  the  element  prestrain  load  vector, 
is  the  element  distributed  load  vector,  and 
is  the  concentrated  gridpoint  load  vector. 

8.  Incremental  Stiffness  Matrix 


The  representation  for  the  frame  element  incorporated  in  the  MAGIC 
System  is  written  to  include  an  incremental  stiffness  matrix.  These  matrices  stem 
from  avoiding  a  complete  linearization  of  the  mathematical  models  for  the  discrete 
elements.  The  formulative  process  is  outlined  below.  Conceptual  examination  of 
this  process  is  deferred  to  the  presentation  of  the  stability  analysis  procedure. 

As  a  first  step,  the  total  strain  induced  at  a  point  is  decomposed  into  a 
contribution  linearly  related  to  displacement  quantities  j-  and  one  which  is  second 
order  in  the  displacement  quantities  ■[  *  N  }  >  i*e’ 

{«T  <  JM€  1  )}+{€N  (  <25> 

Using  this  notation,  the  potential  energy  contribution  which  leads  to  the 
objective  incremental  stiffness  matrix  takes  the  form, 

*o  =  /vl‘<>J  H  {‘N«'}dV  m 
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All  other  energy  terms  associated  with  the  nonlinear  contribution  to  the  total  strain 
are  assumed  to  be  negligible  in  comparison. 

The  knowledge  that  each  term  of  {«(  has  a  linear  dependence  on  the  dis¬ 
placement  functions  and  that  the  dependence  of  each  term  in  /  «N(  )}  is  quadratic 
allows  alternative  expression  of  Equation  26  as 

*c  =  lilj  l  fi  I'gjOhj  (  )dV  <27> 

The  term  Xjj  is  simply  a  multiplicative  constant  and  the  fj,  gj  and  hj  are  dis¬ 
placement  function  forms.  These  displacement  quantities  are  expressible  in  terms 
of  the  assumed  displacement  functions  thereby  accomplishing  the  discretization  of 
the  energy  functional.  Symbolically,  this  expansion  in  terms  of  the  assumed  dis¬ 
placement  function  mode  shapes  can  be  written 


»>JW. 

(28) 

vlB« 

oj{b}j 

(29) 

v  k 

«>]{■}. 

(30) 

The  {b}  matrices  contain  the  independent  variables  of  the  mode  shapes  which 
are  common  to  each  term  of  a  given  element  respresentation  and  the  j  a  j  j, 

(  b|  j,  and  |cj  j  are  the  coefficient  matrices.  The  discretized  potential  energy 
function  which  results  frcu 


from  the  introduction  of  these  assumed  mode  shapes  is  cast 


into  matrix  form 


*crK  IWi-Hj-HJ 


where 


M  'M 


B: 


]T> 

;  x 


M 


■  Lanj 


H 


KJ, dv 
Klj  kb dv 
kb  i< 1  kb  KJj dv 


J.  V' 
1J 


3 

i) 


(31) 


(32) 

(33) 

(34) 
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Now  in  any  given  application  the  vectors  {  a  } 
specified  for  the  contributing  energy  terms.  This 
the  following  items  for  each  term: 

|  b |  j  and  (  c  }  j  must  be 
is  done  by  listing  \ . .  together  with 

Wr  K  J 

tMi 

(35) 

{-WaJ 

1 

(36) 

M.-KJ 

(37) 

The  knowledge  of  these  terms  yields  each  typical  energy  contribution  as 
a  function  of  ‘’he  field  coordinate  displacement  degree s-of-freedom  and  the  sum  of 
the  typical  energy  terms  can  be  carried  out  to  obtain 

*C=T  l  /9j[«]{0}  <38> 

The  matrix  [  N  ]  is  the  element  incremental  stiffness  matrix  referenced 
to  the  field  coordinate  displacement  degrees-of -freedom.  The  previously  defined 
transformation  (Equation  9)  is  introduced  to  obtain  the  element  incremental  stiffness 
matrix  with  reference  to  selected  gridpoint  displacement  degrees-of -freedom,  i.e. 


i  l»q  Jf  N3{»q } 

(39) 

where 

(40) 

The  matrix  [  N  J  is  the  objective  incremental  stiffness  matrix.  It  is 
clear  from  the  foregoing  development  that  the  elements  in  this  matrix  are  functions 
of  the  unknown  displacement  quantities  .  It  follows  that  this  matrix  serves  to 

introduce  the  effects  of  finite  displacements.  The  utilization  of  this  nonlinear  matrix 
is  discussed  as  an  integral  part  of  the  stability  analysis  procedure, 

It  is  recognized  by  the  authors  that  the  foregoing  outline  of  the  development 
of  an  element  incremental  stiffness  matrix  is  lacking  in  clarity.  Matrix  notation  is 
not  well  suited  to  expression  of  nonlinear  relations.  Recourse  to  the  explicit  state¬ 
ment  of  the  incremental  stiffness  matrix  for  the  frame  element  in  Section  4  is  suggest¬ 
ed  for  clarification  of  this  general  symbolic  statement  of  an  element  incremental 
stiffness  matrix. 
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9.  Stress  Matrices 


Having  completed  expression  of  the  total  potential  energy  based  element 
matrices,  it  is  appropriate  to  define  the  element  stress  matrices.  The  element  stress 
matrices  stem  directly  from  the  governing  equations.  The  stress-displacement 
relation  is  obtained  upon  substitution  of  the  strain-displacement  equation  (Equation 
12)  into  the  stress-strain  equation  (Equation  11),  i,e. 

(41) 

Transformation  to  gridpoint  coordinates  and  particularization  to  specific  points  with¬ 
in  the  element  yields 

<42> 

where  the  element  stress  matrices  are  given  by 

[S]  =  [e]  [c]  [r^J  (43, 

{^H  E]{‘i }  <«> 

Stress  resultants  rather  than  point  stresses  are  sought  in  the  thin  shell  and 
slender  prismatic  elements.  Resultants  corresponding  to  deformations  not  consider¬ 
ed  may  be  obtained  directly  from  the  governing  differential  equations  of  equilibrium. 

In  general,  a  rotational  transformation  is  applied  in  order  to  exhibit  stress  values 
with  reference  to  coordinate  axes  which  simplify  interpretation. 

This  completes  statement  of  the  method  employed  in  deriving  the  discrete 
element  representations  incorporated  in  the  MAGIC  System.  The  matrices  of  the 
individual  discrete  elements  are  recorded  in  Sections  4  through  9. 

C.  LINEAR  STRESS  ANALYSIS 

1.  Stiffness  Equation 

The  mathematical  model  for  the  total  structure  is  traditionally  constructed 
by  forming  equilibrium  equations  corresponding  to  the  gridpoint  displacement  degrees- 
of-freedom.  A  more  general  systematic  approach  to  constructing  the  mathematical 
model  for  the  total  structure  is  realized  by  carrying  forward  the  variational  view¬ 
point.  Specifically,  the  energy  functions  for  the  total  structure  can  be  constructed  by 
effecting  a  nonconformable  sum  of  the  individual  element  matrix  energy  forms.  This 
nonconformable  sum,  in  which  common  gridpoint  degrees-of-freedom  are  employed 
for  adjoining  elements,  imposes  continuity  over  the  entire  structure.  Application  of 
the  Euler  equation  (Equation  I)  then  yields  the  objective  governing  equations  for  the 
structure.  This  variational  approach  to  the  assembly  of  elements  to  form  a  total 
structure  representation  is  particularly  attractive  when  generalized  nonphysical  dis¬ 
placements  degrees-of-freedom  such  as  ''wXXy"  are  retained 
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The  element  matrices  for  a  structural  system  are  generated  from  the  in¬ 
put  data  to  the  MAGIC  System  without  regard  to  their  interconnection  as  indicated  in 
Block  2  of  Figure  2.  Since  knowledge  of  these  individual  element  matrices  is  required 
during  subsequent  analysis  phases,  they  are  released  from  the  generation  of  the 
MAGIC  System  as  distinct  submatrices  of  system  level  matrices.  For  example,  the 
linear  strain  energy  stored  in  all  the  elements  is  written  as 


Ul  J  [  KSI  }{ 

where 


A,} 

. K . J 


(45) 


(46) 


(47) 


The  element  column  matrices  are  also  stacked  individually  in  system  level 
vectors.  For  example,  the  external  work  of  all  the  element  loads  applied  to  all 
elements  takes  the  form, 


W  "tAI  ]  {  Pal  } 


(48) 


where 
{  pel 


H(Vfp}i  •{ 


F  +  F  l  - /  F 

*  PJ  2  ’If 


(49) 


Several  additional  system  level  matrices  are  generated  from  the  input  data. 
Firstly,  the  matrix  of  the  gridpoint  loading  at  every  degree-of-freedom  in  all  loading 
conditions  is  provided,  i.e. 
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(50) 


The  input  data  describing  the  interconnection  of  elements  is  processed  to 
obtain  a  system  level  assembly  matrix.  This  assembly  matrix  takes  the  form  of  a 
transformation  matrix  between  all  possible  gridpoint  degrees-of-freedom  |  A  j  } 
and  those  gridpoint  degrees-of-freedom  which  remain  after  interconnection  of  the 


elements  to  form  the  objective  structural 


system  {*.}• 


i.e. 


KMra]{M 


(51) 
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(52) 
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A  further  system  level  matrix  is  generated  to  extract  the  degrees-of- 
freedom  which  actually  exist  from  all  those  associated  with  the  gridpoints  of  the 
assembled  structure.  This  matrix  takes  the  form  of  a  transformation  between  the  com¬ 
plete  set  of  degrees-of-freedom  j  ^  j,  and  the  actual  or  reduced  set  |  ^  ),  ,  i.e. 


where 


KHrJK} 


(53) 


(54) 


Provision  of  the  foregoing  system  level  matrices  enables  execution  of  a 
linear  stress  analysis.  The  first  step  taken  is  to  combine  the  assembly  and  reduction 


transformations  of  Equations  51  and  53  to  obtain 

{A.HFarHAs} 

(55) 

(rarMr.HrJ 

(56) 

This  combined  transformation  is  introduced  into  the  energy  expressions  of 
Equations  45  and  48  to  obtain  the  desired  system  matrices,  i.e. 

Where  <J>p  =  \  [  A  „  j  [  Ks  ]  {  A  s  }  -  [  A  „  J  {  Ps  } 

(57) 

KHrarfK  ][r„  ] 

(58) 

{M  =  [rJT{*>e.Mrr]T{pc} 

(59) 

The  variation  of  this  potential  energy  function  now  yields 
ness  equation  which  takes  the  form 

a  governing  stiff- 

MKHM 

(60) 

This  equation  is  presently  solved  by  inversion.  In  general,  multiple  load  conditions 
exist  and  a  corresponding  multiplicity  of  solutions  is  obtained. 
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2.  Element  Stresses 


As  in  the  case  of  stiffness  and  load  matrices,  the  element  stress  matrices 
are  stacked  individually  in  system  level  matrices,  i.e. 


where 


kHMKl-k} 

Mt-  l }"••"•  M’.-. kfj 

.  K  , 


[SS  ]  = 


•  H 

'  .R 


(61) 

(62) 

(63) 


(64) 


Stresses  are  secondary  variable**  obtained  subsequent  to  the  solution  for 
the  primary  variables  |  As  \  .  Equation  55  enables  direct  expression  of  the  de¬ 

sired  stress  quantities  in  terms  of  the  primary  displacement  variables,  i.e. 


k}=M[rJ{As}-W 


(65) 


3.  Element  Forces 


Element  forces  are  useful  results  in  many  applications.  This  is  particularly 
true  when  the  element  employed  is  a  simulation  of  the  actual  component.  The  back 
substitution  for  element  forces  takes  the  form 


{F„et}  “  (Ks.  ]  [rar]  {as}  -  {peI } 

where 

{  Fnet  }T 

4.  Reactions 


(66) 


The  final  step  is  to  calculate  the  force  balances  and  the  reactions.  These 
are  readily  available  from  tho  element  forces,  i.e., 

T 


D.  STABILITY  ANALYSIS 


An  analysis  procedure  is  incorporated  in  the  MAGIC  System  to  examine  the 
stability  of  flexible  lightweight  structures.  The  structural  stability  phenomena  associ¬ 
ated  with  structures  of  this  type  inevitably  involve  geometrically  nonlinear  behavior, 
while,  in  contrast,  material  behavior  remains  linearly  elastic.  Thus,  it  is  somewhat 
fortuitous  that  geometric  nonlinearities  are  most  readily  incorporated  in  the  subject, 
displacement  methods  of  analysis. 

There  exists  a  hierachy  of  geometric  nonlinearities  which  may  be  incorporated. 
Associated  with  each  level  of  nonlinearity  is  a  degree  of  complexity  and  a  range  of 
applicability.  The  stability  analysis  provided  in  the  MAGIC  System  is  restricted  to 
the  prediction  of  critical  load  values  and  buckling  mode  shapes.  The  prediction  of 
nonlinear  pre-  and  post-buckling  behavior  is  not  attempted. 

The  "classical"  approach  to  buckling  analysis  is  based  upon  the  assumption  that 
the  membrane  force  distribution  induced  in  a  structure  is  known  ab  initio  as  a  linear 
function  of  the  applied  loading.  The  intensify  of  the  given  membrane  force  distribution 
that  causes  the  effective  flexure  stiffness  to  vanish  implies  a  critical  applied  load 
intensity  or  buckling  load. 

The  behavior  of  thin-shell  and  slender  prismatic  structural  components  of 
zero  curvature  is,  within  the  scope  of  linear  mechanics,  naturally  completely  un¬ 
coupled  into  membrane  and  flexure  behavior.  A  similar  uncoupling  of  membrane  and 
flexure  behavior  can  be  obtained  for  components  of  non-zero  curvature  subjected  to 
certain  types  of  boundary  and  applied  load  conditions.  This  uncoupling  of  membrane 
and  flexure  behavior  is  employed  to  advantage  in  the  subject  general  instability 
analysis.  However,  the  nature  of  the  instability  phenomena  associated  with  practical 
built-up  structures  of  complex  configuration  and  applied  loading  transcends  the  scope 
of  the  preceding  classical  buckling  analysis  assumption.  In  general,  membrane  and 
flexure  behavior  can  only  be  uncoupled  within  certain  components  or  zones  of  the 
total  structure.  Linear  coupling,  which  occurs  at  the  junctures,  cannot  be  avoided  and 
is  accounted  for  in  the  subject  general  stability  analysis. 

Geometric  nonlinearities  are  introduced  into  the  analytical  model  for  stability 
analysis  via  the  previously  defined  element  incremental  stiffness  matrix.  In  the 
presence  of  the  incremental  stiffness  matrix,  the  potential  energy  for  a  discrete 
element  takes  the  form 

%  ’i  IsJHKhi-  KJMKHMW  m 

The  element  incremental  stiffness  matrices,  like  the  stiffness  and  load  matrices, 
are  made  available  from  the  matrix  generator  in  the  form 
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The  total  potential  energy  for  a  structural  assembly  of  discrete  elements  can 
be  indicated  symbolically  as 

IAsJKKAsK  [AsJ[Ns]{As}-lAsj{psl  (71) 

in  which  the  definition  of  the  one  new  symbol  introduced,  [  N  j  ,  follows  immediate¬ 
ly  from  the  statement  of  the  linear  stress  analysis,  i.e. 

K]-[r„]T[".,][rJ 

Equation  71  represents  a  geometrically  nonlinear  mathematical  model  suitable 
for  the  prediction  of  certain  types  of  nonlinear  behavior.  The  stability  analysis  in 
the  MAGIC  System  is  directed  toward  the  more  modest  goal  of  predicting  critical 
buckling  loads.  Reference  27  is  recommended  as  a  useful  source  of  information 
regarding  the  prediction  of  nonlinear  pre-  and  post-buckling  behavior. 

The  vanishing  of  the  second  variation  of  the  total  potential  energy  is  invoked  as 
the  buckling  criterion.  Executing  the  second  variation  of  the  potential  energy  of 
Equation  71  obtain, 

K]{SAS}+  KHsaJ-M  (73) 

in  which  the  vector  {  SAs  }  represents  an  arbitary  variation  from  the  displace¬ 
ment  state  /  A  1  . 

lU8  J 

The  computational  utilization  of  Equation  73  for  predicting  critical  loads  is 
based  upon  the  assumption  that  the  incremental  stiffness  matrix  appropriate  to  the 
critical  load  level  pcr  can  be  expressed  as  a  function  of  its  value  at  any  given  load 
level  p  i.e. 


This  is  interpretable  as  assuming  that  the  intensity  of  the  internal  force  state  changes 
linearly  with  changes  in  the  applied  load  without  affecting  the  relative  distribution  of 
the  internal  forces.  Invoking  this  assumption  transforms  Equation  74  into  the  form. 
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(75) 


M  {8 A s}  =  -  (-=-)  [N]  {SAs} 

Clearly,  the  availability  of  a  linear  solution  at  any  load  level  enables  specification  of 
the  matrices  of  this  governing  equation.  The  prediction  of  the  critical  load  level  is 
then  reduced  to  the  solution  of  an  eigenvalue  problem.  The  incorporation  of  this 
approach  in  the  MAGIC  System  provides  a  powerful  tool  for  stability  analyses  of 
general  frame  and  thin  shell  structures. 
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4.  FRAME  ELEMENT 


A.  INTRODUCTION 

A  conventional  frame  discrete  element  is  incorporated  in  the  MAGIC  System. 
This  element,  shown  in  Figure  6,  is  suitable  for  the  idealization  of  structural  compo¬ 
nents  which  are  adequately  characterized  by  "beam  theory."  Having  established 
basic  procedure  and  notation  in  Section  3,  the  mathematical  model  for  the  frame  ele¬ 
ment  is  summarized  in  this  section  in  terms  of  element  matrices.  The  formulation 
is  presented  in  detail  in  Reference  28. 

The  frame  element  is  broadly  applicable  to  space  frame  and  stiffened  shell 
structures.  Connection  eccentricities  can  be  accounted  for  in  shell  stiffener  applica¬ 
tions  of  the  frame  element.  Space  trusses  can  be  accommodated  as  a  special  case  of 
space  frames.  The  truss  specialization  is  particularly  useful  in  combination  with  the 
shear  panel  element  of  Section  5. 

Geometric  specification  of  the  straight  slender  prismatic  frame  element  is 
given,  in  part,  by  the  end  point  coordinates.  A  third  coordinate  point  in  the  positive 
quadrant  of  the  element  axis  system  (Xg,  yg)  is  required  to  specify  the  twist  orienta¬ 
tion. 


The  cross  section  of  the  frame  element  is  assumed  doubly  symmetric  with 
respect  to  the  element  coordinate  axes.  It  is  characterized  by  stiffnesses  I**,  Iyy  and 
Izz  about  the  three  element  axes  together  with  the  cross  sectional  area. 

A  linear  Hooke's  law  is  assumed  to  govern  material  behavior.  Cross  sections 
initially  orthogonal  to  the  element  axis  are  assumed  to  maintain  orthogonality  with  the 
deformed  axis.  It  is  further  assumed  that  deformations  are  sufficiently  small  to  allow 
superposition  of  element  loading. 

Linear  polynominal  axial  and  torsional  displacement  mode  shapes  are  construc¬ 
ted.  A  cubic  polynominal  displacement  mode  shape  is  constructed  for  flexure  in  each 
of  the  two  principal  planes  of  bending.  These  mode  shapes  lead  to  a  total  of  12  unde¬ 
termined  coefficients  for  the  element  which  are  chosen  to  correspond  to  three  trans¬ 
lational  and  three  rotational  displacement  degrees-of-freedom  at  each  end  of  the  ele¬ 
ment.  Description  of  stress  behavior  is  accepted  as  the  definition  of  the  12  forces 
acting  at  the  two  gridpoint  connections. 

Element  matrices  are  provided  for  stiffness,  incremental  stiffness,  stress,  dis¬ 
tributed  loading,  and  axial  thermal  loading.  Certain  of  these  participate  in  the  evalu¬ 
ative  application  to  a  portal  frame  structure  presented  in  this  section. 
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Cross  Section 

Figure  6.  Frame  Element  Representation 
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B.  FORMULATION 

1.  Displacement  Functions 


The  polynomial  mode  shapes  assumed  for  the  displacement  functions  are 

written 


{«  <»}  -  [bo]  {*} 

(76) 

f  I  • 

where 

(u()  }  =  [u(  ),  v(  ),W(  ),$(  )  J 

(77) 

=  [01*02*  ***  *  012  J 

(78) 
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Figure  7.  Displacement  Function  Mode  Shapes  Qb  ( 

and  [bo]  ,  the  matr  ix  of  mode  shapes  is  defined  in  Figure  7.  Elementary  interpo¬ 
lation  theory  is  invoked  to  obtain  a  transformation  to  gridpoint  displacement  degrees- 
of-freedom,  i.e. 


{*}  ■  M{U 


where 


}  -|u,v,w,fl  , .  fl  , ,  Q  ,»u  0 »  v  „ »  w  0»  Q  9  0  o  l(SO) 

1  g  •*  L  gl  gl  gl  °xgl  aygl  °zgl  g2  g2  g2  xg2  yg2  zg2j' 

and  m  is  defined  in  Figure  8. 

Assuming  the  ends  of  the  frame  element  are  positioned  as  shown  in  Figure 
9  relative  to  the  offset  gridpoint,  an  eccentric  connection  transformation  is  provided 


where 


{  Se  )  -luel’vel’wel>  ®xel'  ®yel'  ®zel'ue2’ 


Ve2’  we2’  ®xe2’  ®ye2’  ®ze2j<82> 
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and  [  re]  is  defined  in  Figure  10. 


Since  the  eccentric  connection  is  translational  only,  the  subsequent  trans¬ 
formation  to  system  coordinates  |  xs  |  is  based  upon  the  original  direction  cosines 
of  the  element.  This  transformation  relation  takes  the  form 
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Figure  8.  Displacement  Coordinate  Transformation 


ion 


Figure  9.  Frame  Element  Eccentricity 


where 


The  matrix  (T  |  is  the  direction  cosine  transformation  between  element 
L  gsJ 

and  global  jx  \  coordinate  axes,  i.e. 
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'  9  9  1  9  1  9 

Displacement  Coordinate  Transformation 
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The  final  transformation  builds  in  the  option  to  employ  degrees-of-freedom 
referenced  to  gridpoint  axes.  Transformations  between  the  gridpoint  axes  and  system 
axes  are  known  from  the  input  to  the  MAGIC  System.  For  the  two  gridpoints  of  the 
frame  element  these  take  the  form 


(87) 
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Given  these  axes  transformations,  the  objective  gridpoint  displacement  degrees-of- 
freedom  are  introduced  via  the  relation 


{»s}  ■  [rsqHSq} 

where 


(88) 

(89) 


Q. 


zql*  Uq2*  Vq2*  Wq2’ 


®xq2*  ®yq2’ 


(90) 


In  summary,  displacement  functions  are  assumed  in  terms  of  field  coordi¬ 
nate  displacement  degrees-of-freedom.  These  are  retained  for  the  algebraic  develop¬ 
ment.  Then,  the  foregoing  sequence  of  transformations  is  invoked  to  yield  the  desired 
displacement  degrees-of-freedom.  Collectively,  this  set  of  transformations  is  written 
as 

ffcJK}  (91) 

where 

P>q]  -  [r/9s]  [re]  [rg8]  [rsq]  (92) 

2.  Linearized  Potential  Energy 

The  assumption  of  linear  material  behavior  governed  by 

<r  =  E  ( e  -  t  j)  (93) 

leads  to  a  strain  energy  of  membrane  and  flexure  given  by 

U  =  /  (jE*2  -  E«.*)  dv  (94) 
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The  linearized  strain-displacement  relation  for  the  frame  element  is, 


€  =  u  -  yv  -  zw  (95) 

X  J  XX  XX  v  ' 

The  prescribed  prestrain  6^  is  taken  to  be  constant  over  the  cross  section  and  over  the 
length. 

The  linear  potential  energy  functional  which  arises  in  consequence  of  these 
strain  relations  is 

$pl  =  f  i  (EAu  2  +  El  v  2  +  El  w  2  +  EI  Q  2  (96) 

P1  JK  2  x  z  xx  y  xx  x  x 

-EA€.uv  -  P  V  -  P  W)  dx 

1  X  y  Z 

Note  that  torsion  and  distributed  load  terms  have  been  incorporated  in  the 
above  energy  expression.  The  distributed  loadings  Py  and  Pz  are  assumed  to  be 
linearly  varying  over  the  length  of  the  element,  i.e. 


p<*>=  [i  -  (fr)]pi  +(t)p2 


Substituting  the  displacement  mode  shapes  into  the  linear  potential  energy 
iunctional  and  integrating  over  the  cross-section  obtain, 

where  the  matrices  [k  ],  {  F€}  and  { Fp}  are  given  in  Figures  11,  12,  and  13,  re¬ 
spectively.  In  conformance  with  the  notation  of  Section  3,  the  foregoing  are  the  element 
stiffness  [  K  J  ,  prestrain  load  ,  and  pressure  load  1  Fpj  matrices  referenced 

to  field  coordinate  displacement  degrees -of-freedom. 

It  is  convenient  to  define  a  distinct  prestrain  load  vector  for  strains  induced 
by  thermal  expansions.  Since  flexure  prestrains  are  omitted  in  the  absence  of  know¬ 
ledge  of  the  cro^s  section  geometry,  specialization  of  the  prestrain  j  F*  j  load  to 
thermal  load  j  F€  TJ.  is  accomplished  by  the  relation,  '' 

€i=  a  At  (99) 

3.  Incremental  Stiffness 


The  retention  of  quadratic  displacement  terms  in  the  strain-displacement 
relation  for  the  frame  element  yields  the  strain  contribution, 

. .  1  2  1  2  (100) 
‘N<)  =  2Vx  +  *2  Wx 
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Figure  11.  Stiffness  Matrix  ^KJ 

This  nonlinear  strain  term  leads  to  two  incremental  stiffness  energy  contributions,  i.e. 


where 


(101) 


(102) 


1  2 

—  EAu  w 
2  xx 


dx 


(103) 
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Each  of  these  terms  is  of  the  general  form  constructed  in  the  symbolic  development 
of  Section  3.  The  general  form  is  stated  here  in  the  context  of  the  frame  element  as 


=  r  xfg^.  dx 


where 


>cj = L 

g.  =  [l,x,x2j{b}. 

h.  =  [l,x,  X2J{c}. 


(104) 

(105) 

(106) 

(107) 


The  matrix  form  of  a  typical  contribution  to  the  incremental  stiffness 
energy  now  follows  immediately  in  consequence  of  the  theoretical  development  of 
Section  3,  i.e. 


•„-MWT.WT-M7 


[M 


[cn]  .[Bn] 
’[an] 


(0 

Wj 

Wj 


(108) 


The  three  matrices  [an]  ,  [%]  ,  and  [  ]  are  given  explicit  definition  in  Fig¬ 

ure  14.  Particularization  of  this  general  form  to  the  individual  energy  contributions 
of  Equations  102  and  103  is  given  in  Figure  15  by  specification  of  X  and  the  quantities: 


{a 

}  =N!s,} 

(109) 

(b 

h-  KJ  R}.« 

(110) 

{c 

h-KJl*.}, 

(111) 

The  knowledge  of  these  quantities  enables  the  execution  of  a  nonconformable  sum  of 
the  individual  contributions  to  obtain  the  total  incremental  stiffness  energy  in  terms 
of  the  field  coordinates,  i.e. 

$ 

c 

(112) 

51 


Figure  14.  Incremental  Stiffness  Submatrices  [ANj  «{EN}  and  {Cn} 


The  matrix  [  N  ]  is  the  objective  element  incremental  stiffness  matrix  referenced 
to  field  coordinate  displacement  degrees-of-freedom.  Transformation  to  the  selected 
gridpoint  displacement  degrees-of-freedom  is  accomplished  via  the  previously  de¬ 
rived  transformation; 


{*}  *  [r*]  {\} 


(113) 


Since  the  elements  of  the  incremental  stiffness  are  functions  of  subsets  of  the  {*} 
via  the  {a},  and  {c}j»  they  are  indirectly  functions  of  the  independent  prim¬ 

ary  displacement  variables  {M  .  Thus,  expression  of  the  incremental  stiffness  for 
a  stability  analysis  required  the  availability  of  the  displacement  results  from  a  prior 


analysis. 


The  incremental  stiffness  matrix  for  the  frame  element  has  been  stated 
here  in  accord  with  the  standard  form  outlined  in  Section  3.  Proper  interpretation  is 
very  important  and,  for  this  reason,  a  typical  term  for  the  frame  element  is  resur¬ 
rected  here  and  examined  in  detail. 


52 


4> 


cl 


$c2 


\  =  EA 

Xl  T“ 

{■},-HWT 
Wi =  \ 1,2,3  J  K  05’0r 

H  x =  tt, 


X  =  ^ 

An  O 


0 


1,  2,  3 


^8‘  V  ^10 


Figure  15.  Incremental  Stiffness  Parameters 


Beginning  from  Equation  102,  this  energy  contribution  is  rewritten  as 

<$=/—■  EA  u  v  v'  dx  (1 

cl  J  2  xxx 

to  obtain  the  appearance  of  the  more  general  form  in  which  all  these  functions  are 
different.  This  form  is  rewritten,  in  turn  as 


$ 


cl 


u  V  v' 
xxx 


+  u  v  .v' 
xxx 


+  u  v  v ' )  dx 
xxx 


(115) 


The  over-symbol  requires  definition.  Matrices  are  the  natural  notation  of  multi¬ 
dimensional  linear  algebra.  As  a  consequence  they  have  become  the  language  of  the 
finite  element  analysis  technology.  Recognizing  this,  the  subject  nonlinear  formula¬ 
tion  is  cast  into  matrix  form  to  facilitate  interpretation  and  implementation.  The 
over-symbol  ~  serves  to  identify  variables  which  will  be  imbedded  in  coefficient 
matrices  to  accommodate  matrix  notation.  With  this  explanation,  Equation  115  is 
rewritten  again  as 
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At  this  point  the  assumed  mode  shapes  are  introduced  into  the  vectors  in 
which  ux,  vx,  and  v'  correspond  to  the  f,  g  and  h  of  Equations  105,  106  and  107,  i.e. 


(117) 


Substitution  of  these  mode  shape  quantities  into  the  energy  functional 

yields 
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The  integrated  form  of  this  relation  corresponds  to  the  symbolic  form 
exhibited  in  Equation  108.  The  remaining  step  is  to  bring  in  the  field  coordinate  dis¬ 
placement  degrees-of-freedom.  The  transformation  relation  required  takes  the  form 


(119) 


This  relation  is  equivalent  to  the  symbolic  relations  of  Equations  109,  110  and  111 
taken  collectively.  Introduction  of  this  relation  into  Equation  118  yields  the  objective 
incremental  stiffness  contribution.  These  individual  contributions  are  accumulated 
to  obtain  the  total  as  expressed  in  Equation  112.  Beyond  this  point,  the  symbolic 
statement  of  the  analysis  process  requires  no  further  clarification. 

4.  Stress  Matrix 


Stresses  for  a  frame  discrete  element  are  represented  by  the  end  point 
forces.  These  "stress"  quantities  are  exhibited  with  respect  to  element  axes.  Calcula¬ 
tion  is  based  upon  the  relation 

{  or}  =  [s]{q}  -  {xi}  (120) 

{  a  }  "l  Fxl*  Fyl’  Fzl’  Mxl’  Myl*  M*r  Fx2’  Fy2’  Fz2’  1x2’  My2’  Mz2  J 


where 


Is]  -  [r*J 


(121) 


(122) 


Expression  of  this  stress  matrix  completes  the  specification  of  the  matrices  which 
comprise  the  frame  discrete  element  representation. 
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C.  EVALUATION 


As  an  illustration  of  the  use  of  the  frame  element  in  a  structural  evaluation, 
consider  the  following  problem. 

A  three- member  portal  frame  is  shown  in  Figure  16,  along  with  the  loading, 
pertinent  dimensions  and  material  properties.  The  two  idealizations  used  in  this 
analysis  are  shown  in  Figure  17.  A  comparison  solution  is  given  in  Reference  29 
for  this  portal  frame.  Table  I  presents  the  results  obtained  from  this  analysis  and 
the  reference  solution.  It  should  be  noted  that  the  alternate  finite  element  solution 
neglects  axial  deformation,  thus  producing  a  slightly  stiffer  structure. 


rA 


7 

E  =  10  psi 
y  =  0.30 
A  =  18  in.2 
I0  =  13.5  in.4 

G  =  3.846  x  10  psi 
P  =  550  lb 

Figure  16.  Three  Member  Portal  Frame  Description 


X 


(b)  Six  Element  Portal  Frame 
Figure  17.  Idealizations,  Three  Member  Portal  Frame 


57 


TABLE  I 

COMPARISON  SOLUTIONS  FOR  THREE  MEMBER  PORTAL  FRAME 


Deflection,  pL  (inches) 

Rotation,  6,  x  10-^ 

Node 

Point 

MAGIC 
(3  Elem) 

MAGIC 
(6  Elem) 

Reference 

MAGIC 
(3  Elem) 

MAGIC 
(6  Elem) 

Reference 

2 

0.01143 

-0.7564 

3 

0.02691 

0.02691 

0.02682 

-0.3382 

-0.«382 

-0.3350 

4 

0.02687 

-0.1648 

5 

0.02684 

0.02684 

0.02682 

-0.3376 

-0.3367 

-0.3350 

6 


0.01140 


-0.7544 


5.  QUADRILATERAL  SHEAR  PANEL  ELEMENT 


A.  INTRODUCTION 

A  quadrilateral  shear  panel  is  incorporated  in  the  discrete  element  library  of 
the  MAGIC  System.  This  element,  shown  in  Figure  18,  is  suitable  for  the  represen¬ 
tation  of  thin  membranes  which  carry  load  primarily  by  diagonal  tension.  The  direct 
load  carrying  capacity  of  such  membranes  is  delegated  to  surrounding  axial  force 
members  available  via  the  frame  element  of  Section  4. 

The  general  quadrilateral  shape  of  the  shear  panel  is  defined  by  the  coordinates 
of  the  four  comer  points.  The  geometric  definition  is  completed  by  specification  of 
an  effective  uniform  thickness. 

In  contrast  to  the  usual  approach,  the  principle  of  complementary  energy  is 
employed  to  derive  the  representation  for  the  quadrilateral  shear  panel.  Using  this 
approach,  stress  rather  than  displacement  distributions  are  assumed.  In  particular, 
a  constant  shear  stress  state  is  employed. 

Deformation  behavior  of  the  shear  panel  is  described  by  the  displacements  of 
its  four  corner  gridpoints.  Description  of  stress  behavior  is  accepted  as  the  con¬ 
stant  shear  stress  value. 

The  complete  element  representation  for  the  quadrilateral  shear  panel  is  taken 
to  consist  of  a  stiffness  matrix  and  a  stress  matrix.  These  matrices  are  employed 
in  combination  with  axial  force  members  in  an  evaluative  application  to  a  deep  canti¬ 
levered  beam  in  this  section.  Additional  illustrative  applications  are  included  in 
Section  10. 

B.  FORMULATION 


The  element  representation  for  the  quadrilateral  shear  panel  is  derived  using 
the  principle  of  complementary  energy.  Only  shearing  energy  is  considered.  The 
governing  energy  functional  is  given  by 


(123) 


The  matrices  w  and^0gj-  are  available  from  observation  of  Figure  18,  i.e. 


HT=  F  . 

t  gj  L  2X1’ 


gxl’  Fgyl’  Fgx2’  Fgy2’  Fgx3’  Fgy3’  F»x4’  Fgy4 


(124) 


<».} 


U  V  u  V  .  U  V  U  V  . 

.  gl  gl  g2  g2  g3  g3  g4  g4 


(125) 
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The  shear  stress  function  r  is  chosen  as  the  statically  independent  force 
quantity  and  is  assumed  constant  over  the  element,  i.e. 


r 

xy 


(X,  y)  = 


r  =  Constant. 


(126) 


Recourse  is  made  to  Figure  18  to  obtain  expression  for  the  complete  force 
set  in  terms  of  the  constant  shear  stress.  The  statically  equivalent  corner  grid- 
point  force  set  is  readily  written  as 


h)  ■  k  }  T 


where 


1  <X4  ”X2} 

=L+ — — 


*3  %  t(X4"X2}  +% 

IP  +  2  *  2  ’  2  ’ 


(127) 


(128) 


Substitution  of  Equations  126  and  127  into  the  potential  energy  functional 

yields 


(129) 


At  this  point  an  algebraic  expression  for  the  total  complementary  energy  follows 
by  integration,  i.e. 


*c  -  £  *2  -  T  {rT }  {»,} 

2G 


(130) 


The  variation  of  Equation  129  yields  the  basic  force-deformation  relation  for 
the  shear  panel  in  the  form 


r  -  dik }  k  } 


(131) 


The  introduction  of  Equation  127  now  yields  the  desired  form  of  the  stiff¬ 
ness  matrix  for  the  quadrilateral  shear  panel  referenced  to  the  ( xg,  yg)  coordi¬ 
nate  axis.  This  result  is 


{pg}  ■  [*,]{»,} 


(132) 


61 


where 


KHiH  Mfr. 

Statement  of  the  element  matrices  for  the  quadrilateral  shear  panel  is  com¬ 
pleted  by  introduction  of  coordinate  axis  transformations.  Writing  the  direction 
cosine  transformation  between  element  and  global  axis  as 


} 


(133) 


k)  -  [TJ  k)  <I34) 

leads  immediately  to  the  force  and  displacement  transformations,  i.e. 

k }  -  k  ]  k }  <i35> 


kkkJk) 

The  transformation  matrix  |TgS]  is  given  in  Figure  19.  The  transformation 
to  gridpoint  coordinate  axes  follows  similarly,  except  that  a  distinct  direction 
cosine  transformation  is  associated  with  each  gridpoint,  i.e. 


(136) 


kL  =  kqL  k}j 


(137) 


The  resultant  displacement  and  force  transformations  take  the  form 


[rj  k} 
kk  M  k} 

The  matrix  [rsq]  is  defined  in  Figure  20. 

The  collective  influence  of  the  foregoing  transformation  yields 

M  ■  H  {»,} 

"here  N-  kfkfk]  [rgs][rsq] 


(138) 

(139) 


(140) 


(141) 
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Equation  139  is  the  objective  form  of  the  stiffness  relation  for  the  quadri¬ 
lateral  shear  panel.  Once  the  displacements  have  been  calculated,  the  single  stress 
quantity  r  is  available  from  a  form  of  Equation  130  extended  to  accommodate  the 
transformation  relations,  i.e. 


r  -  [»J  {sq} 

where 

M  -  (n-)  lrrj  [rgs  ]  [rsq  ] 


(142) 


(143) 


This  completes  the  statement  of  the  element  matrix  representation  for  the 
quadrilateral  shear  panel. 


C.  EVALUATION 

As  an  illustration  of  the  use  of  the  quadrilateral  shear  panel  element  in  a 
structural  evaluation  consider  the  following  problem. 

A  cantilever  beam  subjected  to  a  uniform  load  of  2  x  10"  2  kg  per  mm  is  shown 
in  Figure  21  along  with  its  pertinent  dimensions  and  material  properties.  The  three 
idealizations  employed  in  the  finite  element  analysis  are  shown  in  Figure  22.  Axial 
force  members  are  used  along  with  the  subject  shear  panels  for  idealization  of  this 
structure. 

A  solution  to  this  problem  was  obtained  in  Reference  47  utilizing  "equilibrium 
spar  elements"  for  the  web  and  "bar"  elements  for  the  caps.  Solutions  were  also 
obtained  utilizing  spar  elements  with  linear  as  well  as  quadratic  displacement  fields. 

Figure  23  displays  the  tip  and  midpoint  deflections  for  the  three  idealizations 
used  in  this  analysis  along  with  the  displacements  obtained  in  Reference  47.  Figure 
24  shows  the  bending  moment  distribution  for  the  two  element  case  and  Figure  25 
shows  the  moment  dis‘  .ibution  for  the  four  and  eight  element  cases.  It  is  to  be  noted 
that  the  agreement  of  the  bending  moment  distributions  with  the  reference  solutions 
are  excellent. 
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Figure  23.  Tip  and  Center  Deflections  for  Uniformly  Loaded  Cantilever  Beam 


Figure  25.  Bending  Moment  Distribution  for  Four  and 
Eight  Element  Cases 
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6.  TRIANGULAR  CROSS  SECTION  RING  ELEMENT 


A.  INTRODUCTION 

A  triangular  cross  section  ring  element  is  incorporated  in  the  MAGIC  System. 
This  element,  shown  in  Figure  26,  is  suitable  for  the  idealization  of  thick  walled 
axisymmetric  structures  of  arbitrary  shape.  A  detailed  development  of  the  subject 
element  representation  is  presented  in  Reference  30. 

The  ring  element  representation  is  written  with  respect  to  cylindrical  coor¬ 
dinate  axes.  Its  configuration  is  completely  defined  by  specifying  the  radial  and 
axial  coordinates  of  the  three  corner  points.  Anisotropy  is  provided  for  in  the 
physical  and  mechanical  properties  of  the  ring  element.  Orientation  of  orthotropic 
material  axes  in  the  r,  z  plane  is  data  specified. 

Linear  polynomial  functions  are  employed  for  displacement  mode  shapes 
leading  to  constant  element  strain  and  stress  states.  Interelement  continuity  is 
maintained  among  triangular  cross  section  ring  elements  without  explicit  consider¬ 
ation  in  virtue  of  the  straight  edge  displacement  behavior  permitted  by  the  linear 
polynomial  mode  shapes  employed.  The  constant  strain  and  stress  states  within  the 
element  lead  to  a  rc  juirement  for  relatively  fine  idealization  gridworks  when  de¬ 
tailed  stress  behavior  is  desired.  Relatively  coarse  idealization  gridworks  are 
suitable  for  the  prediction  of  stiffness,  displacement  states  and  vibration  character¬ 
istics. 

Distributed  loading  is  assumed  to  exist  against  one  side  of  the  ring  element 
which  provides  for  convenient  consideration  of  pressure  loading.  A  prestrain  load 
vector  is  included  in  the  ring  element  representation  to  accommodate  directly  pre¬ 
strain  and  indirectly  prestress  and  thermal  loading  as  well. 

Deformation  behavior  of  the  triangular  cross  section  ring  element  is  taken  to 
be  described  by  the  six  displacement  degrees-of-freedom  associated  with  the  grid- 
points  which  the  element  connects.  The  stress  behavior  induced  in  the  element 
includes  radial,  circumferential,  axial  and  shear  stress  values. 

Utilization  of  the  triangular  cross  section  ring  element  is  restricted  to  hollow 
structures.  Generally,  simulation  of  a  solid  configuration  can  be  achieved  simply  by 
leaving  a  relatively  small  hollow  cylinder. 

The  discrete  element  technique  was  first  applied  to  the  analysis  of  axisymmet- 
rical  solids  by  Clough  and  Rashid  This  initial  formulation  of  the  triangular 
cross  section  ring  was  extended  by  Wilson  in  Reference  31  to  include  nonaxisymmet- 
ric  as  well  as  axisymmetric  loading.  The  subject  development  follows  Wilson's 
approach,  but  is  restricted  to  the  axisymmetric  case.  The  formulation  is  extended 
beyond  that  of  Wilson,  however,  in  several  ways.  One  of  these  generalizations  is 


71 


that  the  integration  over  the  volume  of  the  ring  is  effected  analytically  under  normal 
circumstances.  Recourse  is  had  to  the  approximate  integration  technique  when  the 
radial  dimension  of  the  ring  is  small  relative  to  the  ring  diameter. 

The  complete  representation  for  the  triangular  cross  section  ring  is  taken 
herein  to  include  matrices  for  stiffness,  pressure  load,  prestrain  load,  thermal  load 
and  stress. 

B.  FORMULATION 

1.  Displacement  Functions 


Linear  polynomial  mode  shapes  are  taken  to  approximate  the  radial  "u" 
and  axial  "w"  displacement  functions  over  the  triangular  cross  section  ring  element. 
With  respect  to  the  global  coordinate  axes  shown  in  Figure  26,  these  mode  shapes 
are  written  as 

{u  ()  }  =  [b()]  {&}  (144) 

where 

{U}T=  [  u  (  ),  w  (  )  J  (145) 


[bo]  =  P’  X’  y’  ’  ’  ] 
L  J  ,  ,  ,1  ,  x , yj 


and  the  are  simply  the  polynomial  coefficients.  These  are  referred  to  as  the 
field  coordinate  displacement  degrees-of- freedom.  Utilization  of  the  linear  assumed 
modes  of  Equation  146  yields  straight-line  edge  displacements  and  assures  satis¬ 
faction  of  interelement  continuity  requirements. 


The  foregoing  assumed  displacement  functions  are  particularized  to 
corner  point  values  and  the  resulting  relationships  inverted  to  yield  a  transformation 
between  the  field  coordinate  displacement  degrees- of- freedom  w  and  gridpoint 
displacement  degrees-of- freedom  .  The  results  are  expressed  collectively 

as 


{*}-  [rssH* } 


where 


kP=[v 


wgl  ’  V  Wg2’  V  Wg3 


The  transformation  matrix  Erjss]  is  defined  in  Figure  27 .  No  further  displacement 
coordinate  transformations  are  required  since  the  .f  8  j  are  suitable  for  assembly 
of  triangular  cross  section  ring  elements.  ® 
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Figure  27.  Displacement  Coordinate  Transformation 
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2.  Potential  Energy 


Linear  elastic  material  behavior  is  assumed.  This  behavior  is  taken  to  be 
governed  by  the  relation 

W-[«]  {{«}-{«,}} 

where 

T 

M  =  [°r  '  aB'  V  °rz  j  (150) 

{‘}T  W  l‘r*  V  Cz-  €rz  j  <151> 

The  matrix  of  orthotropic  elastic  constants  £  E  j  is  specified  in  Figure  28.  The 
vector  |€- j-  is  the  initial  strain  state. 

The  element  potential  energy  is  derived  as  the  sum  of  strain  energy  and 
external  work  contributions.  Invoking  the  stress-strain  relation  of  Equation  148,  the 
strain  energy  is  given  by 

u  =/v  (1/2  M  [E]  {«}  -M  tE]  {*i>) dV  <152> 

In  general,  the  material  property  characterization  is  known  with  reference  to  axes 
orientated  at  an  angle  /m  with  respect  to  the  geometric  axes.  For  this  reason,  it 
is  necessary  to  introduce  a  matrix  for  the  transformation  of  stress  and  strain  states. 
The  desired  transformation  relations  are 

{‘(m>}  -  KJ  {E<g)}  <153> 


These  stress  and  strain  vectors  may  be  interpreted  according  to  Equations  150  and 
151.  The  transformation  ^  J  is  given  in  Figure  29. 


becomes 


Invoking  the  transformation  to  convenient  element  axes,  the  strain  energy 
U  =  /  (1/2  j  [Efe>  ]  {e®}  -  {«.})  dv  (155) 
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Figure  29 .  Stress  and  Strain  Transformation  r€cr 
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where 


M  -  M  [«<“»]  [r€£r] 
{*,}■  ha]"  h,m)] 


The  strain-displacement  relations  appropriate  to  the  axisymmetric  ring 
element  are  given  by 

€r  =  ur  (158) 

€^=  u/r  (159) 

<z  =  wz  (160> 

€  =  u  +  w  (161) 

rz  z  r  ' 

Introducing  the  assumed  displacement,  function  mode  shapes,  the  strains  are  obtained 
in  terms  of  the  field  coordinate  degrees-of -freedom,  i.e. 

{‘<g,"}>  [c <>]{/*} 

where  the  matrix  [c  ]  is  given  in  Figure  30.  This  relation  enables  statement  of 
a  discretized  potential  energy  function  as 

V  l  (1/2i.^J  [c  <  >  ]T  [e(e)  ]  [co  ]{£} 

-  IJ3J  [on]T(r,«})« 

"  J  [  {-  p  sin  a  )u  +  (  p  cos  a  )  w  ]  ~ ■  dr  (163) 

The  last  term  included  in  this  energy  function  is  the  external  work  con¬ 
tribution.  This  arises  in  consequence  of  a  linearly  varying  pressure  distribution 
applied  between  element  gridpoints  1  and  2  as  shown  in  Figure  26.  The  functional 
form  of  this  loading  is 

p  (r,  z)  =  px  +  axr  +  a2z  (164) 


where 


ll  =  ~(riZ2  "Vl)  (P2‘Pl) 
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Figure  30.  Displacement  to  Strain  Transformation  [c()j 


(166) 


Two  algebraic  forms  are  utlimately  given  for  the  pressure  load  vector  to  account  for 
the  special  case  when  is  equal  to 

The  objective  algebraic  form  of  the  total  potential  energy  for  the  triangular 
cross  section  ring  follows  via  integration.  It  is  convenient  to  preface  statement  of  the 
integrated  form  with  the  definition  of  additional  symbolic  notation.  All  integrals 
arising  out  of  Equation  163  are  of  the  general  form 

8..=  f  r*z^dzdr  (167) 

13  4,  r) 

This  symbol  8  jj  is  employed  to  indicate  the  result  of  the  integration.  With  this 
result,  the  integrated  form  of  Equation  163  is  given  by 

®p  - 1/2  l£j  [«]{£} -10]  (5«H 0J  {?p}  <i68> 

The  matrix  [  k]  and  {  Fp  }  ,  given  in  Figures  31  and  32  are  the  objective 
triangular  cross  section  ring  stiffness  and  pressure  load  matrices  referenced  to 
field  coordinate  displacement  degrees -of -freedom.  The  corresponding  prestrain  load 
vector  is  stated  under  the  assumption  of  a  constant  prestrain  over  the  corss  section  as 


where  the  single  new  matrix  [  C  ’  ]  is  given  in  Figure  33.  It  is  convenient  to  have  a 
distinct  load  vector  for  prestrain  due  to  temperature.  The  desired  modification  of 
Equation  166  is 


{«.<m)}  =  AT{a<m)} 


(170) 


where  {  a  m  }  is  the  vector  of  thermal  expansion  coefficients.  It  follows  that  the 
triangular  cross  section  ring  thermal  load  vector,  referenced  to  field  coordinate  dis¬ 
placement  degrees-of-freedom,  is  given  by 


{*„}  =  at  [c']T[t„]  [*<"»]  {„"”>} 


(171) 
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Figure  31.  Stiffness  Matrix  [k] 


"k12  [(Pl+a2  mi2)Sl  +  (ai+a2k12)82] 

"k12  [(Pl  +  a2  ml2)  &2  +  (al  +  a2  kl2)  8  3  ] 

"kl2  [  (P1  +  a2  mi2)mi2  S1  +  (P1  k12  +  mi2  (31  +  2  a2  kl2j)S2  + 
(Pl  +  a2mi2)Sl+(ai  +  a2k12)  ^2 
(pl  +  a2mi2)  82  +  (ai+  a2  ki2)  83 

(P1  +  a2  mi2)  mi2  *1  +  (P1  k12  +  ml2  (ai  +  2  a2  kl2))S2  +  (ai 


}=  2Trl 


[  -[(Pl  +  airl)84  +  a2  85] 

“[(Pl  +  airi)ri  S4  +  a2ri  h] 

-[(Pl+airi)  S5+a2  H\ 


0 

0 

0 


p  (r,  z)  =  Pl  +  V  +  ag  z 


(p2-pl)  Z1 
**  ""riz2“r2zl 


(P2  -  Pl)  rl 
ri  Z2  '  r2  Z1 


Figure  32.  Pressure  Load  Vector  {f  } 
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+  a2  kl2)  k12  83 


82 


0  8  o  o  o  o 

oo  ho  *oi  0  0  0 

o  o  o  o  o  g10 

0  0  sio  0  Sio  0 

Figure  33.  Prestrain  Load  Submatrix  [c] 
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3.  Stress  Matrices 


The  element  stress  matrices  stem  directly  from  the  stress-strain  relation 
of  Equation  149.  The  strains  are  eliminated  from  this  relation  using  Equation  162  to 
obtain  a  set  of  stress-displacement  relations,  i.e. 


{-}=[TC.]T[E(m)]{[Te,][co]{^}-{«r}}  (H2) 


Particularization  of  the  matrix  [C  (  )  J  to  the  centroidal  position  (xc,  yc)  on  the 
element  as  shown  in  Figure  34  yields  the  objective  stress  matrices  for  the  triangular 
cross  section  ring.  Symbolically, 


{cr}=  [?]{*}-  (4  (173) 

where 

[s]=[t<(T]  [is(m)]  [t<0.]  [c(xc,yo>]  (174) 


{4  =  [x„]T  [>>>] 

This  completes  specification  of  the  element  representation  for  the 
triangular  cross  section  ring  element. 


C.  EVALUATION 


As  an  illustration  of  the  use  of  the  triangular  cross  section  ring  element  in  a 
structural  evaluation,  consider  the  following  problem. 

A  thick  walled  circular  disk  in  the  plane  stress  subjected  to  a  radially  varying 
thermal  load  of  the  form  T  =  T0  (1  -  r2)  is  shown  in  Figure  35  along  with  the  load¬ 
ing,  material  properties  and  pertinent  dimensions.  The  three  idealizations  used  in 
the  finite  element  analysis  are  shown  in  Figure  36.  Reference  32  provides  an 
alternative  analytical  solution  for  this  problem  which  is  based  on  the  theory  of 
elasticity.  Figure  37  shows  the  results  for  radial  and  circumferential  stresses  along 
with  radial  displacements  for  the  discrete  element  idealization  shown  in  Figui*e  38  (c). 
Note  that  the  solid  lines  represent  the  alternate  analytical  solution. 
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E  =  1.8  x  107  psi 
v  =  0.30 
T  =  1100  (1  -r2) 


Figure  35.  Thick  Walled  Disk  Subjected  to  Radial  Thermal  Gradient 
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Figure  37 

Stresses  and  Displacements  in  Thermally  Loaded  Disc 


7.  TOROIDAL  THIN  SHELL  RING  ELEMENT 


A.  INTRODUCTION 

A  toroidal  thin  shell  ring  element  is  incorporated  in  the  MAGIC  System.  This 
element,  shown  in  Figure  38,  is  suitable  for  the  idealization  of  axisymmetric  thin 
shells  of  arbitrary  profile.  The  element  configuration  considered  is  that  of  an  arb¬ 
itrary  section  of  revolution  of  a  right  circular  toroidal  shell.  Performance  of  this 
toroidal  ring  element  is  outstanding  relative  to  the  well  known  conic  ring  element. 

The  first  thin  shell  discrete  element  model  put  forward  was  the  singly  curved 
ring  discrete  element  formed  by  a  section  of  revolution  of  a  thin  conical  shell  '  . 

This  element  has  since  been  the  subject  of  numerous  research  investigations  and  re¬ 
ports  35).  The  reasons  for  this  widespread  attention  are  twofold.  Firstly, 

there  exists  a  broad  and  important  class  of  axisymmetric  thin  shell  structures  with 
problematical  axial  variations  which  are  amenable  to  formulation  and  solution  as 
assemblies  of  ring  elements.  Secondly,  behavior  predictions  based  on  the  polygonal 
idealization  afforded  by  the  conic  ring  have  proved,  in  some  cases,  to  be  meaningless. 

Several  papers  haye  attempted  to  lay  down  guidelines  for  avoiding  the  ideali¬ 
zation  pitfalls  (36,  and  for  interpreting  the  predicted  behavior  (3^).  These  papers 
identify  the  primary  sources  of  difficulty  in  using  the  conic  ring  with  the  discontinu¬ 
ities  in  slope  and  stress  which  occur  along  element  circumferential  interface  lines. 
Having  made  this  identification,  it  follows  that  tne  best  response  is  an  element  model 
which  eliminates  the  troublesome  discontinuities. 

Several  discrete  element  models  have  been  reported  which  seek  to  eliminate 
idealization  discontinuities  by  incorporating  curvature  of  the  meridian  in  the  element 
model  These  doubly  curved  elements  have  virtually  eliminated  the  erratic 

stress  predictions  characteristic  of  the  conic  ring.  The  subject  doubly  curved  ring 
element  representation  differs  from  these  primarily  in  the  utilization  of  generalized 
displacement  functions  which  yield  high  precision  stress  predictions.  This  element 
representation  is  developed  in  detail  in  Reference  41. 

The  toroidal  thin  shell  ring  discrete  element  is  formulated  with  respect  to  a 
toroidal  coordinate  system.  In  general,  the  cross  section  profile  of  the  toroidal  seg¬ 
ment  is  circular.  Specialization  to  conic  and  cylindrical  shapes  is  automatically  pro¬ 
vided  for  within  the  MAGIC  System. 

The  geometric  shape  of  the  element  is  specified  by  the  coordinates  and  surface 
orientation  at  its  edge  grid  ring.  The  thickness  of  the  element  is  assumed  constant. 
The  subject  element  is  written  to  accomodate  orthotrooic  materials.  Axes  of  ortho- 
tropy  are  assumed  to  coincide  with  the  principal  axes  of  the  element. 
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The  mathematical  model  for  the  toroidal  ring  embodies  a  coupled  representation 
of  membrane  and  flexure  behavior.  A  state  of  plane  stress  is  assumed  in  formulating 
the  element  representation.  Discretization  is  affected  by  the  construction  of  polynom¬ 
ial  displacement  mode  shapes .  An  osculatory  axisymmetric  polynomial  interpolation 
function  is  taken  to  represent  membrane  displacement  within  the  element.  Trans¬ 
verse  displacement  is  represented  by  a  hyperosculatory  interpolation  function.  Dis¬ 
placement  behavior  is  taken  to  be  described  by  the  ten  displacement  degrees-of-freedom 
which  are  obtained  from  the  polynomial  mode  shapes  at  the  two  grid  rings  connected 
by  the  element.  These  degrees-of-freedom  provide  a  relatively  high  order  of  vari¬ 
ation  in  stress  and  strain  within  the  element.  For  this  reason,  stress  resultants  are 
exibited  at  the  two  boundary  rings  as  well  as  at  the  midpoint  of  the  element. 

The  toroidal  axes  provide  a  suitable  set  of  coordinate  axes  for  assembly  of 
smoothly  connected  toroidal  ring  elements.  If  idealization  discontinuities  are  present 
at  element  junctures,  then  it  is  necessary  to  reference  the  element  representation  to 
a  set  of  global  coordinates.  Global  coordinates  may  be  used  optionally  when  the  tor¬ 
oidal  ring  elements  are  smoothly  interconnected.  The  toroidal  ring  element  is  readily 
specialized  to  yield  end  enclosure  elements.  This  is  a  particularly  useful  feature 
which  was  not  available  in  the  predecessor  co*">.  ring  element. 

The  complete  representation  of  the  toroidal  thin  shell  ring  element  is  taken 
herein  to  include  matrices  for  stiffness,  pressure  load,  thermal  load  and  stress.  The 
toroidal  ring  element  is  somewhat  more  complex  algebraically  than  the  conic  ring 
element.  This  increment  in  complexity  is  given  justification  in  terms  of  improved 
accuracy  with  fewer  elements  in  the  set  of  evaluation  problems  included  in  this  section. 

B.  FORMULATION 

1.  Geometric  Specification 

The  toroidal  shell  parameters  are  obtained  by  reference  to  Figure  38.  The 
basic  coordinate  system  employed  is  toroidal.  This  is  a  right-handed  orthogonal  cur¬ 
vilinear  system.  The  midplane  of  the  shell  is  defined  by  the  (  £  ,  77  )  coordinate  sur¬ 
face.  The  principal  curvatures  of  the  shell  are  aligned  with  the  coordinate  axes. 
Complete  characterization  of  the  system  is  achieved  by  specification  of  the  metric 
parameters  and  the  principal  curvatures. 

The  definition  of  an  element  of  length  (ds)  is 

(d  s  )  2  -  (d  £  )2  +  (d  77  )2  (176) 

where 

d£  -  Ad£  0  <  £  <  a  (  a2-  ax)  (177) 
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(178) 


d  7)  =  B  d/3 

This  leads  immediately  to  the  metric  parameters,  i.  e. 


A  -  1 


3  =  r  + 


[sin  (  a1  +  4)-sinai] 


(  1/a) 


The  principal  curvatures  are  also  found  from  Figure  38, 


p  „  _A_  =  J-  (181) 

£  R£  a 

sin  (  a  +  £/a) 

p  =J_=  - 1 -  (182) 

V  R0  B 

These  expressions  for  the  general  toroidal  configuration  readily  degener¬ 
ate  to  conical  and  cylindrical  ring  cases,  i.  e. 

(a)  Conical  Ring  A  =  1  (183) 

B  =  r,  +  £  cos  a  1  (184) 


'e  =0 


sin  <*. 


Sf£  [<r2'ri)2  +  (Z2- V2] 


(b)  Cylindrical  Ring 


A  =  1 


B  =  r. 


Pt  ~~  ° 

=1/B 

0  <  ^  <  (z  2  -  z  1  ) 
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This  multiplicity  of  parameter  sets  increases  formulative  effort  since 
integrations  must  reflect  the  alternatives;  however,  an  automated  select  feature 
eliminates  any  impact  of  this  multiplicity  in  utilization  of  the  operational  capability. 

The  foregoing  sets  of  parameters,  taken  collectively,  enable  exact  idealiza¬ 
tion  of  cylindrical,  conical,  and  piece  -wise  circular  shells  of  revolution.  More  gen¬ 
eral  shell  profiles  can  be  realistically  approximated  by  combinations  of  these  elements. 

2.  Displacement  Functions 

The  construction  of  admissible  displacement  functions  is  straightforward 
since  the  functions  are  essentially  one-dimensional.  Polynominal  mode  shapes  are 
assumed.  The  membrane  displacement  is  taken  to  be  cubic  in  the  meridianal  arc 
length.  A  quintie  polynomial  is  assumed  for  normal  displacement.  These  assumed 
modes  are  expressed  in  matrix  form  as 

{uo}.[B<,]{/3}  (193) 

where 

{u()}=  [u,  wj  (194) 


r  i  £  £^  £^ 

[b  ( )  J  =  lj  **  ’  ’  ’  3’  4’  5 

L .  ,  ,  .  i.f.t , t  .t  .t  J 


and  the  {*}  are  simply  the  polynomial  coefficients  or,  alternatively,  the  field  coor¬ 
dinate  displacement  degrees -of-freedom.  Transformation  from  the  j)3|  to  gridpoint 
displacement  degrees-of-freedom  is  required  to  enable  proper  interconnection 
with  adjacent  toroidal  thin  shell  ring  elements.  This  transformation  is  effected  by 
imposing  the  following  conditions  on  the  assumed  functions. 


u 

<e> 

=ui 

e=° 

u  (£) 

t =  s 

=  u 

2 

(196) 

ue 

& 

-  u 

£-0  £l 

(£) 

t* s 

=  u 

(197) 

w 

it) 

=  w, 

r°  1 

W  (£) 

taa 

"  w2 

(198) 

we 

it) 

c^o"vei 

W£  it) 

t  s 

-  w 

(199) 

wft 

it) 

w£f 

S 

-  w 

tt  2 

(200) 
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This  set  of  conditions  can  be  expressed  collectively  in  matrix  form  as 


where 


{«}  *  [r  ]  {£} 


{8} 


ul’  u^'  u2  ’  u^2  ’  W1  ’  w£i  ’  "££j,vv2‘  w£2  ’  w££2 


w,  ,  w 


vv 


,  w  9  ,  w  ^  ,  W  , 


(201) 

(202) 


Specific  definition  of  M  is  not  included  since  it  is  the  inversion  of  this 
relation  which  is  desired.  The  inverted  relation  is  written  as 


(203) 


where  the 


Tps  j 11 


is  now  given  explicit  definition  in  Figure  39. 

The  gridpointdegrees-of-freedom  {8}  are  common  to  adjacent  toroidal  thin 
shell  ring  elements  and  are,  therefore,  suitable  for  the  assembly  process.  On  some 
occasions  it  is  convenient  to  use  degrees -of-freedom  referenced  to  a  rectangular 
global  set  of  coordinate  axes.  Moreover,  such  a  system  must  be  employed  if  adjacent 
elements  do  not  interconnect  smoothly.  This  further  transformation  relation  is  given 
by: 


(204) 


where 


{s}“  [rgs]  {a,} 

{sj.  [uj.o.Wj,  0, 

u2,0,w2,0,»(2,0,u(2,  0,  w 
and  the  transformation  [i"1^  ]  is  specified  in  Figure  40. 

The  two  foregoing  transformations  may  be  collected  symbolically  to  obtain 
a  single  transformation  between  the  field  coordinate  land  gridpoint  j  8g  j-  dis¬ 
placement  degrees-of-freedom.  J  J 


(205) 


{0}“  [r/3sKSs}  <206> 

where 

[r/3g)  =  [r/3sHrgs]  l20,> 

This  completes  the  explicit  statement  of  the  displacement  functional  employed  for 
the  toroidal  thin  ^hell  ring  discrete  element. 
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Figure  40.  Displacement  Coordinate  Transformation 


3.  Potential  Energy 

Linear  elastic  material  behavior  is  assumed.  In  accordance  with  this 
assumption,  a  generalized  Hooke's  law  is  employed,  i.e. 


{o  ()}“  [e]  {{«()}-  {«i  (  )}} 

where  „ 

{CT( )}“ 


rEi  =  -  1 

|-E  1/  E  -i 

v  te  e 

( ‘  '"tP  VPi > 

y.  _  E,.  E  ^ 

IP  P'  Pi 

(208) 


(209) 

(210) 
(211) 


The  term 


is  prestrain  state  and  can  be  interpreted  in  accord  with  Equation  210. 


In  virtue  of  the  assumption  of  linear  material  behavior,  the  strain  energy 
can  be  written  as 

»- J[  [*]{«}- M  MM) 

The  next  step  in  proceeding  toward  the  potential  energy  functional  is  to  ex¬ 
press  the  strains  in  terms  of  displacements.  These  equations,  recorded  from  Reference 
42,  are  written 


{«}  -  K,}  +  z  {Af  } 


(213) 


where 


{AnJ  ‘  l  U£  +  X  1  W*  X2U  +  X3W  J 

T 

{Af}  -  [  ,  -X2  J 


(214) 

(215) 


The  quantities  X.  are  defined  as 

(216) 

V'f  X.’  1/Blf-  V'* 

These  arc  given  explicit  definition  by  the  element  configuration  according  to  Equations 
179  through  192. 
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Based  on  these  strain-displacement  relations,  the  total  potential  energy 
functional  is  given  by 

V/0  t’/2  l  i  [  ‘K  ]{Am} 

<•  !/2  [Af  J  [  JK  ]  {Af  }  (21 

-K1W 


-  2  ir  P  w  )  Bd  £ 
z  7 

[  ‘k]  -  2"1  [E] 

{'«}= 2irt  [ E] {*i } 


(218) 

(219) 

(220) 
(221) 


This  completes  the  statement  of  the  potential  energy  functional  for  the  toroi¬ 
dal  thin  shell  ring  element.  The  next  step  in  proceeding  toward  the  objective  element 
representation  is  to  effect  the  discretization  of  the  functional.  Invoking  the  strain- 
displacement  relations  of  Equation  214  and  215  against  the  assumed  mode  shapes  of 
Equation  193  accomplishes  the  discretization  of  the  displacement  functions.  The  results 
may  be  written  symbolically  as 


{Am<>}  =  KojW  <232> 

{A,  (>}-h  (.]{*}  (223) 


The  displacement  to  strain  transformation  matrices 
Figure  41. 


are  defined  in 


The  applied  load  functions  also  require  discretization  by  the  assumption  of 
mode  shapes.  Considering  first  the  pressure  load,  a  linear  variation  is  assumed. 
Translating  this  assumption  into  functional  form  yields 
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o  o 


JL.A  P* 


(224) 


P 

/ 


<£> 


P  ( 0 ) 


(  p(  S)  -  P  (0)  ) 


Compressing  the  notation,  this  expression  is  rewritten 


<£>  - 


J°) 


z 


(225) 


A  similar  linear  form  is  assumed  to  approximate  the  prestrain  load  distri¬ 
bution.  The  corresponding  functional  form  is  given  by 

{«,}  -  {,,<».}  *  z  {«!(»)} 

+  '7  {{ei  ® i  ^ z  ( •§-) {{'C’i**)} ~ }}  (22!i) 

Notational  convenience  is  realized  by  rewriting  this  relation  as 


A  distinct  prestrain  vector  is  provided  in  the  MAGIC  System  for  prestrains 
due  to  temperature.  Specification  of  the  temperature  load  is  accepted  via  for  tempera¬ 
tures,  i.e., 


T 

'2o 


internal  surface  temperature  at  gridpoint  no.  ] 
external  surface  temperature  at  gridpoint  no.  1 
internal  surface  temperature  at  gridpoint  no.  2 
external  surface  temperature  at  gridpoint  no.  2 


The  thermal  prestrain  quantities  follow  immediately  from  this  data, 
quantities  are  defined  as  follows: 

{e.,0)}=  1/2  (tu  +  T10  )  {  <■  } 

h  <10,}=  (*«-  t2„-  Tu- ■>).){«} 

{k.,,0>}=  1/2  (t2.  -  t2o  -th  *tJ{0} 


(228) 

(229) 


(230) 


(231) 
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This  completes  the  definition  of  assumed  functions.  Invoking  these,  obtain 
the  discretized  potential  energy  functional,  i.  e., 

V-C(i/210]MT[ik][d”']  {*} 

.  >/2|^J[d£]T[Jk][d[]{^} 

-l^JK]TK3K0,}-f^Jh]TK]{^n 

-  2ir  (p  <0)+  — p  (1°)  w)  Bd£ 
z  s  z 

Integration  now  yields  the  objective  potential  energy  form  of  the  representa¬ 
tion  for  the  toroidal  thin  shell  ring  element.  The  symbolic  result  is 

-  I#]  {?<}  -  M  {?„}  <233> 

Presentation  cf  these  matrices  is  prefaced  by  definition  of  additional  notation  in  Fig¬ 
ure  42.  Then,  matrices  [k]  ,  |fcJ  and  |Fp  j  ^  a  re  given  explicit  definition 
in  Figures  43,  44,  and  45,  respectively;  The  matrix  m  is  the  element  stiffness 
matrix  referenced  to  field  coordinate  displacement  degrees-of-freedom  {£}  .  The 
matrices  |Fpj-  and  |f€|  are  the  coi’responding  element  pressure  and  prestrain 
loads . 

The  transformation  of  Equation  206  is  introduced  to  reference  thr  element 
matrices  to  degrees-of-freedom  amenable  to  assembly  of  the  elements.  The  result  is 

♦p  ’  1/2  lSd  J  [K]  {*q}  -  l»qj  {Fe}-  [K  j  {Fpl  <224> 

where 

w- 

(vMr if?.} 

{M- 


(235) 

(236) 

(237) 
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Figure  42.  Notation 
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Figure  45.  Pressure  Load  Matrix, 
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This  completes  the  statement  of  the  stiffness  and  applied  load  matrices  for 
the  toroidal  thin  shell  ring  element. 

4.  Stress  Matrices 

An  element  stress  matrix  is  required  to  transform  the  solution  for  the  pri¬ 
mary  displacement  unknowns  to  a  solution  for  the  secondary  stress  resultant  unknowns 
as  well.  Stress  resultants  corresponding  to  deformations  considered  are  available 
directly  from  integrations  of  Equation  208,  i.  e., 


t/3  “  /  dz 
Zff/3 


dz 


(238) 


(239) 


The  calculation  of  shear  stress  resultant  which  is  associated  with  deformation  not 
considered  is  based  upon  equilibrium  requirements,  i.  e., 


Qf  =X2  ♦•>£] 

The  integrated  form  of  the  discretized  stress -displacement  relations  may  be  written 
symbolically  as 


dr. 


dt 


(240) 


{*}  -  [s]  M'M 

where 

{sf- lTf- V  V  V  Q*J  {2i2) 

and  [  s]  is  given  explicit  definition  in  Figure  46.  The  final  form  of  the  clement  stress 
matrix  is  obtained  by  transformation  to  boundary  displacement  degrees -of-freedom,  i.e., 

[s]  ’  [s]  [r«  ]  (2‘l3) 

This  completes  specification  of  the  matrices  which  comprise  the  toroidal 
thin  shell  ring  element  representation. 
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Figure  46.  Stress  Matrix 


C.  EVALUATION 

As  an  illustration  of  the  use  of  the  Toroidal  Thin  Shell  Ring  Element  in  a  struc¬ 
tural  evaluation,  consider  the  following  example  problem. 

A  thin  walled  circular  cylinder,  cantilevered  at  one  end,  is  subjected  to  the 
action  of  bending  moments,  M,  and  shearing  forces,  Q,  both  uniformly  distributed 
along  the  free  edge  of  the  cylinder.  This  cylinder  is  shown  in  Figure  47  along  with 
the  loading,  pertinent  dimensions  and  material  properties. 

Five  finite  element  idealizations  shown  in  Figure  48  were  employed  in  obtain¬ 
ing  results  for  distribution  of  meriodional  moment  in  the  cylinder.  The  results  shown 
in  Figure  49  were  obtained  from  the  16  element  idealization  shown  in  Figure  48  (e). 
Reference  37  provides  an  alternate  analytical  solution  and  it  is  designated  by  the  solid 
line  in  Figure  49. 
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Section  A-A 


Figure  47 .  Cylinder  Subjected  to  End 
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Figure  48.  End  Loaded  Cylinder  Idealizations 
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8.  MALLET  QUADRILATERAL  THIN  SHELL  ELEMENT 


A.  INTRODUCTION 

A  quadrilateral  thin  shell  element  is  incorporated  in  the  discrete  element  library 
of  ihe  MAGIC  System.  This  element,  shown  in  Figure  50,  is  recommended  for  use  as 
the  basic  building  block  for  membranes,  plates,  and  shells.  The  ilclle  triangular  thin 
shell  element  is  a  compatible  companion  element  useful  in  regions  of  irregularity  and 
prominent  double  curvature.  The  Mallett  quadrilateral  thin  shell  element  representa¬ 
tion  is  dcvcloj>ed  in  detail  in  Reference  43. 

The  shape  of  the  general  quadrilateral  element  is  defined  by  the  coordinates  of 
the  four  corner  points.  It  is  a  zero  curvature  element..  The  plane  of  the  element  is 
determined  by  its  first  three  corner  point  coordinates. 

The  subject  element  is  a  thin  shell  element  in  that  both  membrane  and  flexure 
action  are  represented.  Referenced  to  axes  in  the  plane  of  the  element,  the  membrane 
and  flexure  representations  are  uncoupled.  Optional  generation  of  either  or  both  of  the 
representations  is  controlled  by  the  provision  of  associated  effective  thicknesses.  The 
distinct  membrane  and  flexure  effective  thicknesses  are  assumed  constant  over  the 
plane  of  the  element. 


v 

'  o 


Figure  30.  Quadrilateral  Jhin  Shell  Element  Representation 
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Under  normal  circumstances,  four  comer  points  and  four  midside  points  partici¬ 
pate  in  establishing  continuous  connection  of  the  Mallett  quadrilateral  thin  shell  element 
with  adjacent  elements.  Used  in  this  way  input  data  volume  is  reduced  and  accuracy  is 
enhanced.  An  option  is  provided  to  suppress  the  midside  nodes  individually  if  associated 
complexities  arise  in  grid  refinement  or  nonstandard  connections  with  adjacent  elements. 
Invoking  this  suppression  option  causes  linear  variation  to  be  imposed  on  the  specified 
midside  variables. 

The  Mallett  thin  shell  element  is  written  to  accommodate  anisotropy  of  mechani¬ 
cal  and  physical  material  properties.  Orientation  of  material  axes  is  data  specified. 
Temperature  referenced  material  properties,  selected  from  the  materials  library,  are 
assumed  constant  over  the  element. 

A  linear  generalized  Hooke's  Law  is  employed  for  the  equations  of  state.  Three 
options  are  provided;  conventional  plane  stress,  generalized  plane  stress,  and  restric¬ 
ted  plane  strain. 

The  element  formulation  is  discretized  by  the  construction  of  mode  shapes. 
Membrane  displacements  within  the  subject  element  are  approximated  by  quadratic 
polynomials.  Transverse  displacement  is  represented  by  cubic  polynomials.  A  linear 
variation  is  provided  for  midplane  and  gradient  variations  in  thermal  loading.  Other 
element  loadings,  such  as  pressure,  are  assumed  constant  over  the  element.  Deforma¬ 
tion  behavior  of  the  Mallett  quadrilateral  thin  shell  element  is  taken  to  be  described 
by  the  displacement  degrees -of-freedom  associated  with  the  gridpoints  which  it  connects. 

The  variation  in  strain  within  the  element  which  is  permitted  by  the  assumed  dis¬ 
placement  functions,  leads  to  similar  stress  variation.  Advantage  is  taken  of  this  by 
exhibiting  predicted  stress  resultants  at  the  four  corners  as  well  as  at  the  center  of  the 
element.  Inplane  and  normal  direct,  shear,  and  bending  stress  resultants  are  included. 
The  display  of  stress  implies  a  set  of  axes  of  reference.  These  axes  are  data  specified. 


B.  FORMULATION 

1.  Displacement  Functions 

The  displacement  functions  for  the  quadrilateral  thin  shell  element  are 
constructed  with  reference  to  the  oblique  coordinate  axes  (x0»  y<>)  shown  in  Figure  50. 
The  origin  of  this  system  is  taken  at  the  intersection  of  the  diagonals  of  the  quadrila¬ 
teral.  The  orientation  of  these  axes  coincide  with  the  diagonals.  The  xQ  -axis  goes 
through  gridpoint  number  1.  The  yQ  -axis  goes  through  the  first  gridpoint  in  the  counter¬ 
clockwise  direction  which  is  designated  as  number  2. 
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Polynominal  mode  shapes  are  assumed  for  each  of  the  four  zones  shown  in 
Figure  50.  With  respect  to  the  oblique  coordinate  axes  these  mode  shapes  are  written 
as 


0>.  |B  0)  I 

o  L  u  J 

{4} 

j  =  1,  2,  3,  4. 

(244) 

0).  Ib  ®  | 
o  1  v  J 

RJ 

j  =  1,  2,  3,  4. 

(245) 

«>.  |b  0)| 

0  1  w  J 

R} 

j  =  1,  2,  3,  4. 

(246) 

where 

0UT  -  ••••AnieJ  <247> 

te}T  -  Kr  0f2 .  ^n«J  <248> 

and  the  mode  shape  matrices  >  ®}  •  fe,®}  and  {B  «)}  are  given  in  Figures 
51,  52,  and  53,  respectively.  It  isuapparent  from  these  matrices  that  the  mode  shapes 
employed  for  each  zone  are  complete  up  to  the  order  of  truncation. 

Elementary  interpolation  theory  is  invoked  to  obtain  transformation  to 
gridpoint  displacement  degrees- of-freedom,  i.e. 


where 


{*«}- 

[%n  {»«} 

(249) 

tv0]  rj 

(250) 

T  | 

(s„m  j 

\  '  l  “or  “or  “os-  “or  V  “or  "07’  “os- 

Vor  Vo2’  v03'  v04’  v05'  v06’  v07’  Vo8  J 

(251) 

(U 

T  _  1 

lWol*  Wo2’  Wo3’  Wo4’  Woxl*  Wox2’  Wox3’  Wox4* 

Woyl’  woy2’  woy3'  woy4'  wo„5'  wo„6'  wo„7'  wo„8j 

(252) 

and 


are  defined  in  Figures  54  and  55. 


115 


P 

ml 

Pm2  Pm2 

$m4 

0m5  £m6 

Ptm  Pm$  PmB  £mlo  £mll  $ml2  #ml3  $ml4  $ml5  $ml6 

K"'J-  •.  l 

1  , 

x  .  y  . 

2 

x  * 

2 

xy  .  y  » 

•  •III 

.  .  .  .  j 

hi-  [ 

1  . 

x  .  y  . 

• 

2 

xy  »  y  i 

2 

.  X  ,  ,  • 

.  .  .  .  j 

hi-  % 

1  . 

•  x  ,  y  , 

i 

xy  ,  , 

y2  .  x2  , 

.  .  .  .  j 

hi-  1 

1  , 

.  x  ,  y  . 

2 

x  , 

xy  . 

2 

. . 

,  .  .  .  j 

Figure  51.  Displacement  Mode  Shapes  u 

Ptnl 

0  m2  Pm3 

$m4 

Pm5  Ptn6 

$m7  PmB  $m9  £mlO  $mll  $ml2  $ml3  Pml4  Pml5  $ml6 

hi- 

[ 

• 

.  .  .  1  .  X  ,  y  , 

2  2  1 
*  .  xy  ,  y  ,  .  | 

hi-  % 

[ 

*  • 

* 

• 

.  .  .  1  ,  x  .  y  . 

2  2  1 
i  *y  .  y  .  .x 

hi- 

L 

•  • 

• 

.  .  .  1  .  r  ,  y  ■ 

2  2  1 
.  xy  ,  .  y  .  X  J 

hi-  % 

l 

»  t 

. 

. 

>  .  ,  1  ,  x  ,  y  . 

2  2  1 
x  ,  xy  ,  ,  y  ,  j 

Figure  52.  Displacement  Mode  Shapes  v ^ 


Pn 

Pt2 

Pa 

Ph  Pi s  Pk  Pn 

Pn  Pn  Pm  Pin  Pm  Pm  Pm 

Pf\5  #fl6 

[>1  "•  L1  • 

X 

•  y 

X2  ,  2xy  ,  y2  ,  4x3 

,  2  .  2  3 

4x  y  ,  4xy  ,  4y  ,  ,  ,  , 

•  J 

l>"H‘  ■ 

X 

.  y 

.  2xy  .  . 

.23  2 

,4xy  ,  4y  ,  ,  ,  ,  x 

,  3  ,2 

4x  ,4x  y 

lyM* 

X 

.  y 

.  2xy  .  . 

2.2  .3  2 

.  .  .  y  .  4xy  ,  4y  ,  x 

.  3  2 

-lx  ,4x  y 

X 

.  y 

2  „  ,3 

x  ,  2xy  .  ,  ix 

2  2  2  ,3 

4x  y  ,  ,  i  y  .  4xy  ,  4y  , 

■  J 

Figure  53.  Flexure  Displacement  Mode  Shapes  w 


0) 


I —  > 

where 


ft  m9 

ft  ml 

ftmlO 
ft  m2 

ftmll 

ft  m3 

ftml2 

ftm4 

ftml3 

Am  5 

ftml4. 

ftmd 

ftmlo 
ft  ml 

ftmlG 

ftm8 

§9 

*1 

1 

\ 

0 

2 

X1 

0 

0 

0 

0 

&16 

2>2 

1 

0 

y2 

0 

0 

2 

y2 

0 

0 

§n 

1 

X3 

0 

0 

0 

0 

0 

2 

X3 

§12 

»4 

1 

0 

y4 

0 

0 

0 

2 

y4 

0 

[rj- 

§13 

S5 

1 

1 

7xi 

1 

2y2 

1  2 
7xi 

IXiy2 

1  2 
4y2 

0 

0 

§14 

§6 

1 

1 

2  X3 

1 

2  y2 

0 

1  V2 

1  2 
4y2 

0 

1  2 

4  X3 

§15 

§7 

1 

1 

2  X3 

1 

¥y4 

0 

4X3y4 

0 

1  2 
4y4 

1  2 

4  3 

§16 

*8 

1 

1 

2X1 

1 

2y4 

1  2 
4X1 

7xiy4 

0 

1  2 
4y4 

0 

Figure  54.  Membrane  Displacement  Coordinate  Transformation 
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Figure  55.  Flexure  Displacement  Coordinate  IVansforroation 


The  next  step  is  to  introduce  a  transformation  to  gridpoint  degrees-of- 
freedom  referenced  to  the  element  orthogonal  axes  (x  ,  y  ).  The  transformation 
relations  take  the  form 


where 


{*„J  {Sgm} 

KJ  -  [  r.f®3  {>„} 


{  8gm}T  =  i  U£l’  Uc2 . Ug8’  V’  V  '  ’  •’  V  J 


gl’  g2*  1  •  -  g8*  gl’  g2  •  ’  ’  g8 


{  1  T  "l  Wgl’  Wg2’  Wg3’  Wg4’  V*’  Wgx2’  Wgx3’  W 


w 


gyl*  Wgy2’  ”gy3’  "n5’  ”n6*  "n7 


w 


w 


w 


w 


r-wnJ 


gx4’ 


(253) 

(254) 

(255) 

(256) 


and  the  transformation  relations 
and  57. 


are  given  in  Figures  56 


At  this  point  the  optional  transformations  for  the  suppression  of  midpoint 
displacement  degrees-of-freedom  are  introduced.  This  feature  provides  flexibility 
in  idealization  and  facilitates  eccentric  connection  of  elements.  The  transformations 
take  the  form, 


{\»}  {  k } 

{*gE}  =[C°] 


The  degrees-of-freedom  {8’gm}  and  {§’  f }  may  be  interpreted  accord¬ 
ing  to  Equations  255  and  256.  The  suppression  trans  formations  [rsup<m>]  and 
[  rgUp(^]are  defined  in  Figures  58  and  59.  The  effect  of  these  transformations  is  to 
build  in  linear  variations  in  place  of  the  degrees-of-freedom  suppressed. 
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Figure  57.  Flexure  Displacement  Coordinate  Transformation 
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Figure  ">9.  flexure  Displacement  Coordinate  Transformation 
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At  this  point  a  transformation  is  defined  to  establish  a  vectorial  sign  con¬ 
vention  for  the  rotational  degrees-of-freedom  associated  with  the  four  corner  grid- 
points.  In  addition,  the  midside  rotations,  if  not  previously  suppressed,  are  assigned 
to  be  vectorially  positive  from  the  comer  point  with  the  smaller  gridpoint  number 
toward  the  comer  point  with  the  larger  gridpoint  number.  This  transformation  is 
written  symbolically  as 


sgn 


gf 


(259) 


where 


{sgf}  =  lv  V  V  V  8gxl.  egx2.  egx3.  0gx4. 

8gyl'  ®gy2’  8gy3’  Sgy4’  ®n5”  ®n6’  ®n T  ®n8  J 
and  the  transformation  matrix  [  Tggn'^  ]  is  exhibited  in  Figure  60. 


(260) 


A  second  set  of  optional  transformations  is  introduced  to  enable  eccentric 
connection  of  the  quadrilateral  thin  shell  element  to  a  surface  which  is  a  distance  e 
above  the  element.  This  transformation  takes  the  form 


{»„>■ 


(261) 


The  degrees-of- freedom 


{» 


gm. 


H 

rix  [r, 


l} and 

Equation  255,  and  the  transformation  matrix 
that  utilization  of  this  transformation  requires  the  presuppression  of  midpoint  dis¬ 
placement  degrees-of-freedom. 


may  be  interpreted  according  to 
e  ]  is  shown  in  Figure  61.  Note 


Global  or  ’'system'*  displacement  degrees-of-freedom  are  obtained  by  the 
introduction  of  further  transformations  of  the  form: 


(262) 


where 


{8;Hrg/>]  W 


(263) 


Figure  60.  Flexure  Displacement  Coordinate  Transformation 
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Figure  61.  Eccentric  Connection  Transformation.  [1^ 
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(U  l  Usl*  Vsl*  Wsxl*  ^syl*  ^szl*  Us2’ Vs2*  Ws2’  ^sx2’  ^sy2’  ^sz2* 
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(264) 


The  transformation  matrices  .  [Vm,l  and  £  rgs(f)  ]  are  given  in  Figures 
62  and  63.  Assembly  of  the  thin  shell  element  can  be  referenced  to  these  system  dis¬ 
placement  degrees-of-freedom;  however,  it  is  convenient  in  many  cases  to  employ 
special  gridpoint  coordinate  axes.  Accordingly,  a  final  transformation  to  gridpoint 
displacement  degrees-of-freedom  is  provided,  i.e. 


W-  [rj  W 

where 
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(265) 


(266) 
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Figure  62.  Membrane  Coordinate  Transformation 
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The  matrix  [  F  Sqj  is  made  up  of  the  individual  gridpoint  axes  direction 
cosine  transformations  from  the  relations 

{Xs}=[Tsq],  {xq}j  ,267' 

J 

positioned  along  the  major  diagonal  as  shown  in  Figure  64. 

The  foregoing  transformations  may  be  collected  symbolically  to  obtain  a 
single  transformation  between  the  field  coordinate  displacement  degrees-of-freedom 
j  0  ■  ".nd  the  final  gridpoint  displacement  degrees-of-freedom  jq  j  .  The  results 
are  a»  ivilows:  1  J 


iV  -  [Vm,]M 

(268) 

(M  -  [Vf)]  (4 

(269) 

where 

rr  <m)'  n 

[I>q,m,]  =  [r>o(m>]  [Vm>]  [rsup,m)] 

!  F(m)  1 

V-i....  <: 

Lie  J 

:r„.m 

(270) 


[vf,i  ■  V]  [roB(f>]  [rsup<£)]  rcf)]  iv'}  trJ  (27i> 


This  completes  the  explicit  statement  of  the  displacement  functions  employed  for  the 
Mallett  quadrilateral  thin  shell  element. 

.2.  Potential  Energy 

The  strain  energy  density  for  a  thin  shell  element  of  zero  curvature  is 

defined  as 


where 


dV  m  J |d<  J  {o’) 

(272) 

HT  -  Iv  V  vJ 

(273) 

'T' 

£ 

(274) 
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Given: 


Known: 

Plane  Stress  cr  =  0  > 

Z  1 

>  by  definition 

Plane  Strain  «  =  0  J 

2 

Observe: 


Figure  64,  Plano  Stress/Strain  Option 
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Linear  elastic  material  behavior  is  assumed  to  take  place  from  an  initial 
state  of  strain  jgjJ  to  a  final  state  of  stress  |  <7  j.  and  strain  j  e  j  , 

{o-*”'}7  ,  [e<“>]  {{<<”>}-  {««}}  (275) 

The  matrix  of  elastic  constants  [  E  J  is  given  explicit  expression  for  the 
special  case  of  orthotropv  in  Figure  64.  The  superscript  m  indicates  the  coordinate 
axes  of  reference. 

Substitution  of  the  assumed  constituitive  relation  into  the  strain  energy 
density  definition  yields,  after  integration,  an  expression  for  the  strain  energy  of 
elastic  deformation  in  terms  of  the  strains. 

dU  =  i  [£'m>j  [><“>]  le<m>]  (276) 


If  the  material  axes  (m)  are  orientated  at  an  angle  with  respect  to  the 
chosen  element  geometric  axes  (g)  a  transformation  must  be  introduced, 

{«(m)}  =  K<r]  {«(g>}  <277) 

{cr(g)}  =  [tco.]  (278) 

The  transformation  [  T<<r  j  is  defined  in  Figure  65.  Transforming  the 
axes  of  reference  of  the  strain  energy  aensity  and  the  constituitive  equation  obtain 


dU  =  A  l«»J  [■»]  {.»}  -  l««J  {»«} 

and 

{<rfe,}  =  [e«>]  {«®} -{«,<">} 

where 

[■*]  -  M T  [»(m|]  [t„] 

{«t<g>}=[^]T  tE<m>]  w6>} 


(279) 

(280) 

(281) 

(282) 
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[V] . 


€  2  2 
x  +  cos  y  ,  +  sin  /,  +  sinycos  y 

€  2  2 
y  +  sin  y  ,  +  cos  y  ,  -sinycos/ 

<  2  2 
xy  -2  sinycos  y,  +2  sinycos  y,  +  cos  y-  sin  y 


Figure  65.  Strain  Transformation  ^f0-J 


The  well  known  strain-displacement  relations  for  a  thin  shell  element  of 
zero  curvature  can  be  written  as  a  sum  of  membrane  {  e  j  and  flexure  {  K  }  contri¬ 
butions. 


{«<*>}  =  {e»>}  ♦  z  {Kfe>} 


It  is  convenient  to  separate  the  membrane  strain  into  linear  and  nonlinear  parts. 


Explicit  definition  of  the  strain  contributions  in  terms  of  the  displacement  is  given  by 
T 

^  t _ \  * 


{.«}  -I 

u  , 

V  ,  u  +  V  I 

l  u  J  L 

X 

y  y  x  J 

kte)}  =11 


12  12  I 

—  w  ,  —  w  ,  W  W  I 
2  x  2  y  x  y  J 


{K®}  =l-w 


,  -W 

XX  yy 


It  is  convenient  to  carry  forward  the  separation  of  membrane  and  flexure 
strains  into  the  strain  energy  expression  prior  to  introducing  these  strain-displace¬ 
ment  relations.  In  so  doing  it  is  assumed  that  |  j  is  a  linear  function  of  the  z 
coordinate. 


{«»}”  {‘mi}"  2  {«fl} 
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The  resulting  expression  for  the  strain  energy  is  written  as  a  sum  of  mem¬ 
brane  4*^,  ,  flexure  $£  ,  and  coupling  &c  contributions.  Including  an  external  work 
term  <J>p  as  well,  the  set  of  four  energy  contributions  is  written  as, 

*m  -  /  (!k<E)J  f^8’]  k(8)}  leu®J  {?mi<B)})dA  (289) 

"A 

*t  -  f  (I3  Lk<s>J  [£<g)]  {K<g)} -SlK<g)J  KJ) «  i29o> 

w  A 


■  j[  (*  l»  «]  [Efe>]  {..«})  dA 

"A 

*P  ■  /  k»(g))dA 


This  set  of  energy  functions  is  employed  as  the  point  of  departure  in  deriving 
the  companion  triangular  thin  shell  element  representation  in  Section  9. 

In  order  to  realize  the  algebraic  simplification  afforded  by  oblique  coordi¬ 
nates  axes,  it  is  necessary  to  transform  the  displacement  functions  of  Equations  285, 

286  and  287  before  substituting  into  the  strain  energy.  Given  the  transformation  relation 


,  cos  on  c 
’  ”  sin  a  I  Xg 


OJ  L  ’  sinaJl/gJ 

and  using  the  chain  rule  for  partial  differentiation,  the  following  transformation  rela¬ 
tions  are  derived. 


f  =  f 
x  x 
g  O 

(  =  -°2°2.  i  +  4—  I 

y  sin  a  x  sin  a  y 

^  /v  n  r 


f  f 

XX  =  XX 

g  g  0  0 
2 

„  cos  a 


x  y„ 
g  e 


.2  XX 

sm  a  00 

COStl  ^ 

sin  a  x  x 
0  0 


2  y  y 
sin  a  o  0 


+  -7 —  f 

sin  a  xv 
00 


2  cos  a 
2 

sm  a 


x  y„ 
oJo 


134 


Invoking  these  transformation  relations,  285,  obtain  the  strains  expressed  in 
terms  of  displacement  functions  defined  with  reference  to  the  element  oblique  coordinate 
system: 


c 

5? 

n 

I-u3  {  Amu  } 

(299) 

k“}  ■ 

[T  1  { A  } 

L  wJ  v  mwJ 

(300) 

|K'g)  }  = 

[TJ  {Aft,  } 

(301) 

where 

{Amu}T 

II 

(302) 

T 

(a  } 

l  mwJ 

112  12 
=  —  w  ,—  w  ,  w  w 

12  x  2  y  xyJ 

(303) 

kfW } T 

=  -  w  ,  -w  ,  -2w 

L  xx  yy  xy  J 

(304) 

The  matrices  [  Tu  j  and  [  Tw  ]  are  given  in  Figure  66. 

Introducing  these  strain-displacement  relations  obtain  the  energy  functionals 
in  terms  of  displacements  referenced  to  the  oblique  coordinate  system. 

n 

to  |  1-1 

L  AmuJ  [imk]  {AmJ  -  [AmuJ  {imf})  dxd. 

(305) 

*,  -/  (i 

*/xy 

lAJ  [’fij]  {Af }  -  [Af  J  {lfJ  )  dxdy 

(306) 

*.  ■  /  (I* 

•'xy 

muJ  [‘cl  {Amw}) 

(307) 

%  -  J  (px 

*xy 

where 

w  sin  a  j  dx  dy 

(308) 

[  *mk  ]  tsina 

k]V]  k] 

(309) 
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-w  w 
yy  2  y 


-2  w  ,  w  w 
xy  x  y 


12  1  2 
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Figure  66.  Displacement  Function  Transformations 
T 

••  _  *  ~  / _ \  « 


{lm(}=  tai-a  [tJ  {?m«>} 

[■J  =4^K]t  K]  [\] 

o  T 

fT  \  _  t3  sing  r  -I  r-  (g)\ 


[■J 

t  sin  a 

r. 

12 

1 

t3  sin  a 

h 

tw 

12 

l1 

h) 

=  t  sin  a 

[T, 

(311) 

(312) 

(313) 


Equations  303,  304,  305  and  306  are  the  desired  form  of  the  energy  functional. 
It  should  be  noted  that,  in  expressing  the  nonlinear  coupling  energy,  the  second  order 
terms  in  the  prestrain  have  been  assumed  small  relative  to  the  corresponding  first 
order  terms  in  the  total  potential  energy  functional. 
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The  next  step  in  constructing  the  element  representation  is  to  effect  the 
discretization  by  introducing  the  previously  derived  mode  shapes  into  Equations  300, 
303,  and  304.  This  results  in  the  relations. 


{Amul  [DmU)]  0U  <314> 

<3i5> 


The  matrices  [>“]  and  [Df0)]  are  presented  in  Figures  67  and  68, 
respectively.  The  vector  j  Amw}  is  a  quadratic  function  of  the  coordinates  |)9  f  } 
and  symbolic  representation  is  not  attempted  at  this  point. 


Algebraic  statement  of  the  membrane  energy  contribution  of  Equation  305 
is  considered  first.  Examination  of  the  component  relations  of  Equation  314  leads  to 
identification  of  a  typical  form  for  each  element  of  the  vector  (A  mu}  > 


KJ/’  -  M  kJ/>  {<■„}/' 


(316) 


where 

{d}T  =  [l.x.y  J 


(317) 


For  example,  focusing  on  the  first  zone  (k  =  1),  the  first  element  (£  =  1)  is  given  by 


(318) 


Explicit  statement  of  the  [  cm  J  and  {  Ctm|  matrices  for  each  of  the  four  zones  is 
given  in  Figures  69  through  72. 
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1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 


x 


Q.  1.  0,  2x,  y,  0,  0,  0, 


°»  o,  1,  0,  x,  2y,  0,  0, 


,  0,  1,  0,  2x,  y,  0,  o,  0 


,  0,  0,  1,  0,  x,  2y,  o,  0 


v  . 
x 


°»  0,  0,  y,  0,  o,  2x, 

0,  0,  l,  0,  x,  2y,  0,  0, 


0,  1,  0,  0,  y,  0,  0,  2x 


.  0,  0,  1,  0,  x,  2y,  0,  0 


0,  1,  0,  0,  y,  o,  0,  2x, 

°»  0,  1,  0,  x,  0,  2y,  0, 


0,  1,  0,  0,  y,  0,  0,  2x 


,  0,  0,  1,  0,  x,  0,  2y,  0 


y  L 


o,  1,  0,  2x,  y,  0,  0,  0, 


0,  0,  1,  0,  x,  0,  2y,  0, 


0  >  •  1 »  0 , .  2x )  y ,  0 ,  0 ,  0 


,  0,  0,  1,  0,  x,  0,  2y,  0 


Figure  67.  Membrane  Displacement  Derivative  Matrices 
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l 


Figure  70.  Zone  2  Membrane  Stiffness  Parameters  Figure  71.  Zone  3  Membrane  Stiffness  Parameters 
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The  general  form  identified  with  the  elements  of  the  vector  {Amu}  leads 
naturally  to  a  general  form  for  the  associated  energy  contributions.  Firstly,  the  mem¬ 
brane  energy  is  expressed  in  indicia!  notation. 


k  -  i  i  x  (oj: 

k  =  l  j - 1  i  =  l 


(k) 

ij 


(319) 


where 


(k)  ,  ,  (k)  (k) 

-f  (j  (mk)ij  (\u).  (Amu  ) 

^  •'zone  k  1  ^ 


-(In,e).(Amu),<k))dXdy 


(320' 


The  general  contributing  energy  form  now  follows  directly  by  introducing 
the  general  .form  for  the  elements  defined  in  Equation  316,  i.e. 

(k)  .  (k)  (k)  (k)  (k)  (k) 

(*»)..  M,  h»J.  K]  U.  {U 


ij 


ij 


j 


(k)  (k)  (k) 

),KJ1  f°nj  ,  {°<} 


where 


Kj"  -  X0nek  ld 


dy 


dx  dy 


(321) 


(322) 

(323) 


Presentation  of  these  matrices  is  prefaced  by  definition  of  notation  in 
Figure  73.  Then,  [c^J  and  {c^jare  given  in  Figure  74.  The  knowledge  of  these  mat¬ 
rices  together  with  that  of  the  J" omJ  and  {®m}  matrices  specified  in  Figures  69  through 

72  enables  explicit  algebraic  expression  of  each  of  the  ($m  These  individual  energy 
contributions  are  summed  to  obtain  the  objective  algebraic  expression  for  the  total  mem¬ 
brane  potential  energy,  i.e. 


I’m  -rKj  [*„]  OU  -  KJ  {?<} 


(324) 
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Figure  74.  Stiffness  Submatrices 


As  disclosed  by  the  notation  employed,  [  Km]  and  {  Fg}  are  the  membrane  stiffness 
and  prestrain  matrices  referenced  to  field  coordinate  displacement  degrees-of-free- 
dom  w.  .  Explicit  statement  of  these  matrices  is  redundant  since  they  are  simply 
the  assembled  results  of  explicitly  specified  contributions. 

"It  is  convenient  to  have  separate  load  vectors  for  prestrains  due  to  temper¬ 
ature.  The  desired  modification  is  available  immediately  from  Equation  310. 

{l€}  =ATmtsina[Tu]  [Te<r]  [E,m)]  {a(m)}  (325) 

The  objective  algebraic  statement  of  the  flexural  potential  energy  follows 
in  analogy  with  the  development  for  the  membrane  potential  energy.  Examination  of 
the  component  relations  of  Equation  315  leads  to  identification  of  a  typical  form  for 
each  element  of  the  vector  |  A  f  J  ,  i.e. 

r  i(k)  i  i  r  i(k)  r  i(k) 

tv,;  -  M  hb  Kb  (326) 


For  example,  focusing  on  the  first  zone  (k  =  1),  the  first  element  1)  is  given  by 
(1) 


(V 


i'JL 


w 


(1) 


XX 


x  ,  yj 

-2, 

,  -24, 

- 

#f7  ' 

,  -8 

^f8 

to  to 

(327) 


Explicit  statement  of  the  f  cf  J  and  f  }  matriee  for  each  of  the  four 
zones  is  given  in  Figures  75  through  78. 
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Figure  75.  Zone  1  Flexure  Stiffness  Parameters 


(Af).  ‘ 

- 1 

a 

t 

- 1 

5° 

to 

00 

L, - 

[•.] 

It 
►  = 

^f!4’  ^fl5’  ^fl6  J 

(Ar)2  ; 

i 

[  2,  8,  24  J 

\ 

M 

T 

r- 

[**’  *»’  ^fioj 

> 

CO 

°f 

= 

4,  16,  16  J 

\ 

T 

>  = 

_  ^f5’  ^fl6’  J 

Figure  76.  Zone  2  Flexure  Stiffness  Parameters 
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Figure  77.  Zone  3  Flexure  Stiffness  Parameters 
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Figure  78.  Zone  4  Flexure  Stiffness  Parameters 
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The  general  form  identified  with  the  elements  of  the  vector  1“  £ |  leads 
to  a  general  form  for  the  associated  energy  contributions.  As  before,  the  energy 
function  is  expressed  in  indicial  notation  as 


where 


4  3  3 

111 

=  i  ;  =  i  i  -  i 


(328) 


(329) 


The  general  contributing  energy  form  now  follows  directly  by  introducing  the 
general  form  for  the  elements  |A  '  j,  defined  in  Equation  326,  i.e. 


r  (U  l-,JW  f°f J 

'  ij  i 

,  ,  {k)  r  ,  (k) 

(hi).  lafJ.  fcfJ. 


(330) 


Particularization  of  this  general  form  to  the  individual  [  C£  J  and  {(*  f  J 
and  summation  yields  the  objective  algebraic  expression  for  the  total  flexural  poten¬ 
tial  energy,  i.e. 


(331) 


As  disclosed  by  the  notation  employed,  Kf  ]  and  {  }  are  the  flexure  stiffness 

and  prestrain  matrices  referenced  to  field  coordinate  displacement  degrees-of- free¬ 
dom  1*1  .  Explicit  statement  of  these  matrices  is  omitted  since  they  are  simply  the 
assembled  results  of  explicitly  specified  contributions. 


As  in  the  case  of  the  membrane  prestrain  load  vector,  the  flexure  prestrain 
load  vector  is  particularized  to  thermal  loading.  The  desired  modification  follows 
immediately  from  Equation  309,  i.e. 


{‘J  =ATt 


t3  sing 
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Work  equivalent  gridpoint  forces  are  provided  for  the  case  of  a  transverse 
load  uniformly  distributed  over  the  quadrilateral  thin  shell  element.  The  external  work 
of  this  loading  is  defined  by 


/  p*  w 
►  XV 


sin  a  dx  dy 


The  introduction  of  the  assumed  displacement  modes  into  this  expression  yields 
4 

^p  =  Z  /*  Pz  sin  a  [/?fj  {Bw}(k)dxdy 
k  =  1  'zone  k 


Substitution  from  Equation  246  and  integration  then  yields  an  algebraic  expression  for 
the  external  work,  i.e. 


-  l0(J  00 


p  L^fj  l  pJ 

The  matrix  ,  referred  to  as  the  pressure  load  matrix  referenced  to  field 

coordinate  displacement  degrees-of-freedom,  is  given  in  Figure  79. 

3.  Stress  Matrices 

The  stress  resultants  for  a  thin  shell  of  zero  curvature  are  defined  in  the 
notation  of  Figure  80  as  follows: 
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Figure  79.  Pressure  Load  Vector, 
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Figure  80.  Stress  Resultants 


It  was  tacitly  assumed  in  defining  the  stress  resultants  that  nor '  Inear  mem¬ 
brane  flexure  coupling  contributions  to  the  stress  resultants  are  small  relative  to  first 
order  terms.  This  assumption  is  carried  forward  in  writing  the  stress  resultants  in 


terms  of  the  strains 

{n*>}.  t[E»] 

K]Ku<0,}-‘{ 

?  *  t  to-  te)} 

mi  '  l  mo  J 

(340) 
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[■"](*«") -4 

[G2]  [E'g>]  {dyf K<S>' 

}  (342) 

where  [«]-  [l:V> 

)■[“]-[;  t:l 

(343) 

The  stress  resultants  are  expressed  with  reference  to  displacement  functions 
defined  in  the  oblique  coordinate  system  of  the  element  by  substituting  from  Equations 
299,  300,  and  301. 
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Introducing  the  displacement  mode  shapes  assumed  over  the  four  zones  of 
the  element,  the  stress  resultants  can  be  written  collectively  as 
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C.  EVALUATION 


1.  Membrane  Stress  Analysis 

The  first  illustration  which  uses  the  quadrilateral  thin  shell  element  in  a 
structural  evaluation  will  be  the  following.  Consider  a  thin  square  isotropic 
plate  loaded  with  a  self  equilibrating  parabolic  membrane  load  as  shown  in  Figure  81. 
The  material  properties  and  pertinent  geometric  data  are  also  shown  in  the  figure. 

The  idealizations  used  for  the  finite  element  analyses  are  shown  in  Figure  82. 
Three  different  grid  sizes  were  employed  in  this  evaluation.  One  element,  four  element, 
and  16  element  solutions  were  obtained  in  order  to  evaluate  convergence  characteristics. 
Due  to  conditions  of  symmetry  it  should  be  noted  that  only  one  quadrant  of  the  plate  was 
analyzed.  For  the  finite  element  idealizations  employed  in  this  evaluation,  the  midside 
nodes  which  were  loaded  by  the  parabolic  load  were  suppressed.  This  suppression 
invokes  a  linear  edge  displacement  under  the  load. 

The  results  obtained  from  this  set  of  convergence  studies  are  presented  in 
Figures  83  and  84.  Figure  83  is  a  plot  of  the  membrane  displacement,  uq,  at  the  middle 
of  the  plate's  edge  versus  degrees-of-freedom  employed  in  the  analyses.  The  reference 
solution  (Reference  44)  is  designated  by  the  solid  line.  Figure  84  presents  a  curve  of 
the  membrane  displacement,  ux,  and  stress  resultant,  Nx,  versus  the  edge  span  of  the 
plate  for  the  idealization  shown  in  Figure  82  (16  element  solution). 


Figure  81.  Parabolically  Loaded  Membrane 
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Figure  83  Membrane  Displacement  At  Midside  vs  Degrees  of  Freedom 
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Figure  84  Membrane  Displacement  And  Stress  Behavior 
Versus  Plate  Edge  Span 
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Again  it  should  be  noted  that  the  reference  solution  is  designated  by  the 
solid  lines  and  no  discernible  difference  between  solutions  can  be  detected. 

2.  Membrane  Gridwork  Influences 

The  second  illustration  which  utilizes  the  quadrilateral  thin  shell  element 
in  a  structural  evaluation  will  be  the  following.  Again,  consider  a  thin 
square  isotropic  plate  loaded  with  a  self  equilibrating  parabolic  membrane  load  as 
shown  in  Figure  81 .  This  illustration  will  involve  the  effect  that  the  shape  of  the 
elements  used  in  the  structural  idealization  has  on  the  determination  of  the  center 
edge  displacement  of  the  membrane. 

The  six  idealizations  used  for  the  shape  study  are  shown  in  Figure  85. 

It  should  be  noted  that  due  to  symmetry  only  one  quadrant  of  the  plate  was  analyzed. 

The  midside  nodes  which  were  loaded  by  the  parabolic  membrane  load  were  suppressed 
in  this  solution. 

The  results  obtained  from  the  subject  shape  studies  are  shown  in  Figure  86. 
These  solutions  indicate  that  the  displacement  values,  uq,  obtained  for  the  middle  of 
the  plate's  edge,  are  fairly  insensitive  to  the  shape  of  the  element  for  this  class  of 
problem. 


3.  Plate  Stress  Analysis 

The  third  illustration  which  utilizes  the  quadrilateral  thin  shell  element  in 
a  structural  evaluation  v/ill  be  the  following.  A  simply  supported  isotropic  square  plate 
with  a  uniform  normal  pressure  load  of  one  psi  is  shown  in  Figure  87  along  with  its 
material  properties  and  pertinent  dimensions. 

The  idealizations  used  for  the  finite  element  analyses  are  shown  in  Figure 
82.  Note  that  no  node  points  are  suppressed  in  this  analysis.  Three  different  grid 
sizes  were  employed  in  this  evaluation.  One  element,  four  element  and  16  element  solu¬ 
tions  were  obtained  in  order  to  evaluate  convergence  characteristics.  Due  to  condi¬ 
tions  of  symmetry  only  one  quadrant  of  the  plate  was  analyzed. 

Figure  88  is  a  plot  of  the  transverse  displacement  at  the  center  of  the 
plate  versus  degrees-of-freedom  employed  in  the  analyses.  The  reference  solution 
(Reference  45)  is  designated  by  the  solid  line.  Figure  89  presents  a  curve  of  the 
transverse  displacement,  wx,  and  bending  moment,  Mx,  versus  the  center  span  of  the 
plate  for  the  idealization  shown  in  Figure  82  (c)  (16  element  solution).  Again,  it  should  be 
noted  that  the  reference  solution  is  designated  by  the  solid  lines  and  no  discernible 
difference  between  solutions  can  be  detected. 

4.  Plate  Gridwork  Influences 

The  fourth  illustration  which  utilizes  the  quadrilateral  thin  shell  element  in 
a  structural  evaluation  will  be  the  following.  A  simply  supported  isotropic  plate  with  a 
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Figure  87.  Simply  Supported  Square  Plate  with 
Uniform  Normal  Load 
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Transverse  Displacement,  w_  x  10  (in) 


Legend: 

O  Magic  Quadrilateral 

____  Reference  Alternate  Solution 


Figure  88  Convergence  of  Quadrilateral  Element 

Plate,  Square,  Isotropic.  Simple  Support,  Uniform  Pressure 


a  uniform  normal  pressure  load  of  one  psi  is  shown  in  Figure  87  along  with  its 
material  properties  and  pertinent  dimensions. 

This  illustration  will  involve  the  effect  that  the  shape  of  the  elements  used 
in  the  structural  idealization  has  on  the  determination  of  the  maximum  displacement 
of  the  plate. 

The  six  idealizations  used  for  the  shape  study  are  shown  in  Figure  85. 

Due  to  conditions  of  symmetry  only  one  quadrant  of  the  plate  was  analyzed. 

The  results  obtained  from  the  subject  shape  studies  are  shown  in  Figure  90. 
These  solutions  indicate  that  the  determination  of  the  plate’s  center  transverse 
displacement  is  fairly  insensitive  to  the  shape  of  the  element  for  this  class  of  problem. 


161 


162 


0.  IIKU.K  TRIANGULAR  THIN  SHELL  ELEMENT 


A.  INTRODUCTION 

A  triangular  thin  shell  element  Ik  incorporated  in  the  discrete  element  library 
<if  the  MAGIC  System.  This  element,  illustrated  in  Figure  91 .  is  recommended  for 
use  as  the  basic  building  block  for  most  doubly  curved  shells.  Additionally,  it  is 
useful  in  combination  with  the  Mallett  tiuadrilater.il  thin  shell  element  for  dealing 
with  irregular  geometries  of  all  membrane,  plate,  and  shell  structures.  The  Helle 
triangular  thin  shell  element  representation  is  develop'd  in  detail  in  Reference  40. 


Figure  !>1.  Triangular  Thin  Shell  Element  Representation 


The  shape  of  the  general  triangular  element  is  defined  by  the  coordinates  of  the 
three  corner  points.  It  is  a  zero  curvature  element.  The  plane  of  the  element  is 
determined  by  the  three  corner  point  coordinates. 

The  subject  element  is  a  thin  shell  element  in  that  both  membrane  and  flexure 
action  are  represented.  Referenced  to  axes  in  the  plane  of  the  element,  the  membrane 
and  flexure  representations  are  uncoupled.  Optional  generation  of  either  or  both  of 
the  representations  is  controlled  by  the  provision  of  associated  effective  thicknesses. 
The  distinct  membrane  and  flexure  effective  thicknesses  are  assumed  constant  over 
the  plane  of  the  element. 

Under  normal  circumstances,  three  corner  points  and  three  midside  points 
participate  in  establishing  continuous  connection  of  the  Helle  triangular  thin  shell 
element  with  adjacent  elements.  Used  in  this  way  input  data  volume  is  reduced  and 
accuracy  is  enhanced.  An  option  is  provided  to  suppress  the  midside  nodes  indi¬ 
vidually  if  associated  complexities  arise  in  grid  refinement  or  nonstandard  con¬ 
nections  with  adjacent  elements.  Invoking  this  suppression  option  causes  linear 
variation  to  be  imposed  on  the  specified  midside  variables. 

The  Helle  thin  shell  element  is  written  to  accommodate  anisotropy  of  mechani¬ 
cal  and  physical  material  properties.  Orientation  of  material  axes  is  data  specified. 
Temperature  referenced  material  properties,  selected  from  the  Materials  Library, 
are  assumed  constant  over  the  element. 

A  linear  generalized  Hooke’s  Law  is  employed  for  the  equations  of  state. 

Three  options  are  provided;  namely,  conventional  plane  stress,  corrected  plane 
stress,  and  restricted  plane  strain. 

The  element  formulation  is  discretized  by  the  construction  of  mode  shapes. 
Membrane  displacements  within  the  subject  element  are  approximated  by  quadratic 
polynomials.  Transverse  displacement  is  represented  by  cubic  polynomials.  A 
linear  variation  is  provided  for  midplane  and  gradient  variations  in  thermal  loading. 
Other  element  loadings  such  as  pressure  are  assumed  constant  over  the  element. 


Deformation  behavior  of  the  Helle  triangular  thin  shell  element  is  taken  to  be 
described  by  the  displacement  degrees -of -freedom  associated  with  the  gridpoints 
which  it  connects. 

The  linear  variation  in  strain  within  the  element  which  is  permitted  by  the 
assumed  displacement  functions  leads  to  similar  stress  variation.  Advantage  is 
taken  of  this  by  exhibiting  predicted  stress  resultants  at  the  three  corners  as  well 
as  at  the  center  of  the  element.  Inplane  and  normal;  direct,  shear,  and  bending 
stress  resultants  are  included.  The  display  of  stresses  implies  a  set  of  axes  of 
reference.  These  axes  are  data  specified. 
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B.  FORMULATION 


1.  Displacement  Functions 

The  displacement  functions  for  the  triangular  thin  shell  element  are  con¬ 
structed  with  reference  to  the  coordinate  system  (xg,  yg)  shown  in  Figure  91.  The 
origin  of  this  system  is  located  at  the  centroid  of  the  triangle.  The  orientation  of 
the  xg  axis  is  defined  by  corner  gridpoint  number  1.  The  yg  axis  of  this  right- 
handed  coordinate  system  is  taken  counterclockwise  from  the  Xg  axis  in  the  plane 
of  the  element. 


Polynomial  mode  shapes  are  employed  to  represent  the  displacement 
functions  over  the  element.  These  mode  shapes,  for  the  membrane  displacements, 
take  the  form 


where 


“  =  lBu  J{*m} 

*-h  JOU 


{am  }  |_0ml  1  $m2 . 0ml2  j 

The  mode  shape  matrices  [  Bu  J  and  [bv]  are  exhibited  in  Figure  92. 
be  noted  that  both  of  these  represent  complete  quadratic  polynominals. 


(353) 

(354) 

(355) 
It  should 


The  transverse  displacement  function  is  approximated  by  distinct  polynomials 
over  each  of  the  three  zones  of  the  triangular  element  identified  in  Figure  91.  The 
basic  cubic  polynomials  may  be  written  symbolically  as 

w(j)  =  Ia  J  {r(J)}  j=i-2-3-  (356) 

where  the  mode  shape  matrix  [l3wJis  given  in  Figure  92  and  the  ^y^^are  simply 
the  undetermined  coefficients. 


Interzone  continuity  requirements  impose  interdependencies  among  these 
undetermined  coefficients,  yielding 


{rU)}  -  M  {*,} 

The  resulting  admissible  displacement  functions  are  given  by 

(•TkJ  M  W 

where 

{^f  }  =  l^fl  ’  ^f2  ’  '  '  *  ’  ^fl2  J 

The  transformation  to  independent  field  coordinate  displacement  degrees- 

of-freedom  for  transverse  displacement  \PA  are  exhibited  in  Figures  93  through  95. 
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Figure  92.  Displacement  Mode  Shape  Matrices 
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Figure  93.  Flexure  Displacement  Interzone  Transformation 
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Figure  94.  Flexure  Displacement  Interzone  Transformation 


Elementary  interpolation  theory  is  invoked  to  obtain  transformations  from  field 
coordinate  to  gridpoint  j  S  j  displacement  degrees -of -freedom,  i.e. 
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It  is  convenient  to  define  the  transformation  matrices 
in  terms  of  submatrices  by  writing 


-1 


(364) 


(365) 


These  submatrices  are  now  given  explicit  definition  in  Figures  96  and  97. 

At  this  point,  optional  transformations  are  introduced  to  enable  suppression 
of  the  midside  displacement  degrees -of-freedom.  This  feature  provides  flexibility 
in  idealization  and  facilitates  consideration  of  eccentric  connections.  The  transfor¬ 
mations  take  the  form 

{8gm}-  fW"’]  {8'gm}  <366) 

Km)  =  kp®]  {8’gJ  <36,> 

The  degrees -of-freedom  {$'gm}  and  {S'gf}  may  be  interpreted  according  to  Equa¬ 
tions  362  and  363.  The  suppression  transformations  f  T  and  [  T  (f)]  are 

given  in  Figures  98  and  99.  *■  sup  ■*  *•  sup  •* 


A  transformation  related  only  to  flexure  is  defined  to  establish  a  vectorial  sign 
convention  for  the  rotational  degrees -of-freedom  associated  with  the  three  corner 
gridpoints.  In  addition,  the  midside  rotations  are  signed  so  as  to  be  compatible  with 
adjoining  elements.  Specifically,  they  are  assigned  vectoriaily  positive  in  the  direc¬ 
tion  from  the  corner  gridpoint  with  the  smaller  gridpoint  number  toward  the  corner 
gridpolnt  with  the  larger  number.  This  transformation  relation  is  written  symbol¬ 
ically  as 
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Figure  97.  Flexure  Displacement  Cooriinatt  Transformation 
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Figure  98.  Flexure  Displacement  Coordinate  Transformation 


The  transformation  matrix  [  *s  shown  in  Figure  100. 

The  transformation  introduced  next  is  designed  to  enable  eccentric  connection 
of  the  triangular  thin  shell  element  to  a  surface  which  is  a  distance  ez  above  the  ele¬ 
ment.  This  transformation  effects  a  coupling  of  the  flexure  degrees -of -freedom  to 

the  membrane  displacements  and  is  written  as 

•  ^ 

{*'«■}-  [M-M]  '  f  ^  ■  <370> 

[{*"«}  J 


173 


174 


wgl  gxl  gyl  Wg2  gx2  gy2  Wg3  gx3  gy3  n4  ®n5  ®nG 


gl 


gxl 


gyi 


g2 


gx2 


gy2 


g3 


gx3 


gy3 


n4 


+1 


-1 


,  “I 


+1  , 


1  , 


,  +1 


"I  , 


,  ±1  , 


n5 


,  ±1  , 


n6 


,  ± 


Figure  100.  Flexure  Displacement  Coordinate  Transformation 

The  degrees -of -freedom  {&"gm}  and  {8"gf  }  may  be  interpreted  according  to 
Equations  362  and  369  with  the  understanding  that  these  quantities  now  refer  to  the 
eccentric  gridpoints.  The  eccentric  connection  transformation  [r.]  is  given  in 
Figure  101.  Note  that  utilization  of  this  eccentric  connection  feature  requires  the 
presuppression  of  the  midside  gridpoint  displacement  degrees -of-freedom. 

Global  or  "system"  displacement  degrees -of-freedom  are  introduced  via  a 
further  set  of  transformation  relations  of  the  form 


{sv-HVI  K> 
{8"gf}  =  M  W 


(3 
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The  matrices  [V  *]  and  £  rgs(f)]  which  accomplish  the  transformation  to 
system  displacements  are  given  in  Figures  102  and  103.  These  system  displace¬ 
ment  degrees -of-freedom  can  be  employed  to  assemble  discrete  elements;  however, 
in  many  applications,  it  is  convenient  to  employ  special  gridpoint  coordinate  axes. 
Accordingly,  a  final  transformation  relation  is  provided  to  reference  the  displacement 
degrees -of-freedom  to  gridpoint  coordinate  axes,  i.e. 

W-MW  <3,4> 

where 

KJMvvv  V  V  V  (375) 

uq2'  V  w<i2’  ®qy2' 

Uq3’  Vq3’  Wq3’  ^qx3’  ^qy3*  ^qz3’ 

V  V  V  ®qn4'  °  •  °  ' 

V  V  V  V  °  •  °  ' 

V  V  V  ®qn6'  °  ’  °  j 

The  matrix  [  Tgg  ]  is  made  up  of  the  individual  gridpoint  axes  direction 
cosine  transformations  from  the  relation 


K}*  [Tsq]  k}j 


positioned  along  the  major  diagonal  as  shown  in  Figure  104. 

The  foregoing  transformations  may  be  collected  symbolically  to  obtain  a  single 
transformation  between  the  field  coordinate  displacement  degrees -of-freedom  j  j 
and  the  final  gridpoint  displacement  degrees -of-freedom  {Si.  The  results  are  as 
follows:  ’  ' 


{AJ  “  [^q(m)]  {»q} 

kl  -  M  {»,} 
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Figure  103.  Flexure  Displacement  Coordinate  Transformation  T 


Figure  104.  Displacement  Coordinate  Transformation 


] 


This  completes  the  explicit  statement  of  the  displacement  functions  constructed  for 
the  Helle  triangular  thin  shell  element. 


2.  Potential  Energy 


The  total  potential  energy  functional  appropriate  to  thin  shells  of  zero  curva¬ 
ture  was  stated  in  Equations  289  through  292  in  terms  of  four  contributions,  i.e. 


-  j[(i kS)J  kg)]  k^Hkl  {‘Jg)})  « 


(381) 


-  /  (t (fr)  [k<b)J  [E<gl]{K<g)}-(fe)lK<s)j  l\(g)})  dA  <382) 
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*c  -  X  (‘Kte>J  [E<E)]  k(g)}) dA  (383) 

A 

--  J  (pzw(g))  dA  (384) 

Accepting  the  foregoing  statement  of  the  total  potential  energy  function  as 
the  point  of  departure,  the  first  step  is  executed  primarily  for  notational  convenience. 
The  strains  are  written  in  terms  of  the  displacement  functions  via  the  relations 


fe®}T-(A  }TfT  1T=|u,u,v,v  I 
l  U  J  l  muj  L  Uj  L  X  y  x  yj 


{ew<E’}T=  (Amw}T  =li'WV  WxWyj 

{K<e,}MAfw}T5l-'V  -V2wxyJ 


1,  0,  0 
0,  o,  1 
0,  c,  1 
0,  1,  0 


(385) 


(386) 

(387) 


The  potential  energy  functional  contributions  for  the  triangular  thin  shell  element 
are  now  written  in  analogy  with  those  for  the  quadrilateral  thin  shell  element  as 


<t>  =  f  (±|A  Ifl  }-Ia  I  (i  })dA 

m  J  \2L  muJLmkJl  mu  J  L  muj  l  mcJ  / 
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®f  /*  (t  l^fwj  [jfk]{^fw}  "  l^fwj  {!f«r} )  dA 
JA 

*c  =  J(lAnJ[ln][*J)dA 

A 

<X>  =  /  fp  w|dA 

P  JA  \  z  I 

where 

fImkj't[Tu]T[T<<y]T[E<m’]K<r]  [TJ 

l^mc}  =  l[Tu]T  [T<c ]  ^  [E<m)] 

W  4KJtIe<m,]KJ 


(388) 
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(393) 

(394) 
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The  next  step  in  constructing  the  element  representation  is  to  effect  the  dis¬ 
cretization  by  introduction  of  the  previously  derived  displacement  mode  shapes  into 
Equations  385,  386,  and  387.  This  substitution  yields 

{Amu}  -  [Dm]  {*0  <397> 

{AtHDf]  {y(k)}  <398» 

The  matrices  ^  j  and  W  are  presented  in  Figures  105  and  106.  The  vector 
{ Amw}  *s  a  quadratic  function  of  the  coordinates  |  y  j  and  symbolic  representation 
is  not  attempted  at  this  point. 

Algebraic  statement  of  the  membrane  energy  contribution  of  Equation  388 
is  considered  first.  Examination  of  the  component  relations  of  Equation  397  leads  to 
identification  of  a  typical  form  for  each  element  of  the  vector  j  Amjl ,  i.e. 


where 


(Amu)^=  M  KnJ, 

Mt  =  I1*  y>J 


(399) 


(400) 


For  example,  focusing  upon  the  first  term  {/=  1)  obtain, 

(Amu)/=\=ll'X-y  J  n  1  ^ 


1, 


t  2, 


,  1 


m2 
ft  m4 

$m5j 


(401) 


Explicit  statement  of  all  of  these  [cm  J  and  [<*m  j  matrices  is  given  in  Figure  107. 

The  general  form  identified  with  the  elements  of  the  vector  {  A  mu  }  leads 
naturally  to  a  general  form  for  the  associated  energy  contributions.  Firstly,  the  mem¬ 
brane  energy  functional  is  rewritten  as 
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Figure  106.  Flexure  Displacement  Derivatives  Matrix  |D^ 
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Figure  108.  Definition  of  Integral  Notation, 
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The  general  contributing  energy  form  now  follows  directly  by  introducing 
the  general  form  for  the  elements  defined  in  Equation  401,  i.e. 

V^m)ij  _  3  (Imk  )ij  [am]i  [°mj  i  [°mk]  [Cmjj  {Qm}j 

(404 

-('nub  |°mji[0mji{°m«} 


where 


[cmk  ]  =  -/*{d}Ld  J  dXdy 

[°m«]'  /{d}dxdy 


Presentation  of  these  matrices  is  prefaced  by  definitions  of  notation  in  Figure  108. 
Then,  [cmk]  and  [cm<]  are  given  in  Figure  109.  The  knowledge  of  these  matrices, 
together  with  the  statements  of  the  [cmk]  and  {am}  matrices,  enables  explicit 
expression  of  each  of  the  (^m)  jj.  These  individual  contributions  are  summed  to 
obtain  the  objective  algebraic  expression  for  the  total  membrane  potential  energy,  i.e. 


■  7  KJ  pm]  {4}-  m  {?m«} 


The  matrix  J  is  the  element  membrane  stiffness  matrix  and  the  vector  {?»«} 
is  the  element  membrane  prestrain  load  representation.  As  disclosed  by  the 
notation,  both  of  these  are  referenced  to  the  field  coordinate  displacement  degrees - 
of -freedom.  Explicit  statement  of  these  matrices  is  considered  redundant  since 
they  are  simply  the  assembled  results  of  known  contributions. 

It  is  convenient  to  define  a  special  prestrain  vector  for  thermal  loading. 
This  is  easily  accomplished  by  rewriting  Equation  393  to  read 


{\.«}-AT«‘  [Tu]T  [T«r]  T  [E<m)]  [«<m)] 


The  objective  algebraic  statement  of  the  element  flexural  potential  energy 
follows  in  analogy  with  the  development  for  the  element  membrane  potential  energy. 
Examination  of  the  component  relations  of  Equation  398  leads  to  identification  of  a 
typical  form  for  each  element  of  the  vector  f  j-,  i.e. 
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Figure  109 .  Definition  of  Notation, 
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For  example,  focusing  on  the  first  zone  (k  =  l),  the  first  element  (X-  1)  is  given  by 


(410) 


Explicit  statement  of  all  these  and  [a*]  matrices  is  given  in  Figure  110.  It 

should  be  noted  that  these  matrices  are  common  to  the  three  zones  of  the  element. 

Having  defined  the  foregoing  quantities,  the  flexure  energy  is  accountably 
rewritten,  i.e. 


*£  -  Z  l  l  (#f)„ 

k  =  l  i  =  l  j  =  1 '  x/y 


(k) 


(411) 


where 


{♦«)„*  - /  HU/’ (*,).* 


zone  k 


(412) 


The  general  contributing  energy  form  now  follows  directly  by  introducing 
the  general  form  for  the  elements  of  ^  A  defined  in  Equation  410,  i.e. 

hLhr 
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Figure  110.  Definition  of  Notation,  [ C  Jand  {ajj  } 

Particularization  of  this  general  form  to  the  individual  [c{  J*  {  «c} 
and  ^  matrices  and  summation  yields  an  algebraic  expression  for 
the  element  flexural  potential  energy  associated  with  the  kth  zone,  i.e. 

M(k,{M<k) 


(414) 


Substitution  from  Equation  357  brings  in  the  field  coordinate  displacement  degrees - 
of-freedom 

(415) 

A  final  summation  over  the  three  zones  yields  the  flexure  stiffness  [  Kf  ]  and  pre¬ 
strain  load  {Ff<}  matrices  referenced  to  field  coordinate  displacement  degrees -of- 
freedom.  This  result  is  written 
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As  in  the  case  of  the  membrane  prestrain  load  vector,  the  flexural  prestrr.in 
load  vector  is  particularized  for  thermal  loading.  The  desired  modification  follows 
immediately  from  Equation  395,  i.e. 


{U  “AT£l2  M  [E<m)]  {“ 


(417) 


The  final  element  matrix  which  arises  from  the  potential  energy  is  the  pres¬ 
sure  load  matrix.  The  external  work  term  of  Equation  389  is  expanded  in  terms  of  the 
displacement  mode  shape  for  the  zone  to  yield 


w=p,  /  K’J  K}  Akt 

H 


(418) 


Integrating  this  expression  obtain 

W  .  ly«/>j  {*F/>} 


(419) 


wht.>  is  given  in  Figure  111.  Summation  over  the  three  zones  is  accomplished 

with  the  introduction  of  Equation  357.  The  results  an  expression  for  the  external 
work  which  contains  the  pressure  load  vector  /  Fp\  referenced  to  field  coordinate  dis 
placement  degrees -of-freedom  ^  J 

w  '  KJR}  (42°) 
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Figure  111.  Pressure  Load  Vector 


3.  Stress  Matrices 


The  stress  resultants  for  a  thin  shell  of  zero  curvature  are  defined,  in  the 
notation  of  Figure  80,  as  follows: 


N  =  /  <r  dz  N  =  /  r  dz  (421) 

y  J  y  xy  J  xy 

z  z 


189 
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It  was  tacitly  assumed  in  defining  the  stress  resultants  that  nonlinear 
membrane-flexure  coupling  contributions  to  the  stress  resultants  are  small  rela¬ 
tive  to  first  order  terms.  This  assumption  is  carried  forward  in  writing  the  stress 
resultants  in  terms  of  the  strains 

K®} ■= 1  [e<b)]  k®}  -  *  k®}  <424» 

K®}  -  4  [E®]  k®}  -4  K®}  <425> 

k®} -  4  M  k®]  {Jr k®}} 4N  kfe,]{iHK<B)}} 
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H  ■ 

0,  0,  ll 

_0, 1,  oj 

The  stress  resultants  are  expressed  in  terms  of  the  displacement  func¬ 
tions  by  substitution  from  Equations  383  and  384. 
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Introducing  the  displacement  mode  shapes  assumed  over  the  three  zones  of 
the  triangular  thin  shell  element,  the  stress  resultants  can  be  written  collectively  as 
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This  completes  the  statement  of  the  matrices  which  comprise  the  Helle 
triangular  thin  shell  element. 


C.  EVALUATION 


1.  Membrane  Stress  Analysis 

The  first  illustration  which  utilizes  the  triangular  thin  shell  element  in  a 
structural  evaluation  will  be  the  following.  A  thin  square  isotropic  plate  loaded  with 
a  self  equilibrating  parabolic  membrane  load  is  shown  in  Figure  112  along  with  its 
material  properties  and  pertinent  geometric  data. 


t  =  0.1  in. 

E  =  30  x  106  lb/in.2 


v  =  0.3 


Figure  112.  Parabolical ly  Loaded  Membrane 


The  idealizations  used  for  the  finite  element  analyses  are  shown  in  Figure 
113.  Three  different  grid  sizes  were  employed  in  this  evaluation.  Two  element, 
eight  element  and  32  element  solutions  were  obtained  in  order  to  evaluate  convergence 
characteristics.  Due  to  conditions  of  symmetry  only  one  quadrant  of  the  plate  was 
analyzed.  For  the  finite  element  idealizations  employed  in  this  evaluation,  the  mid¬ 
side  nodes  which  were  loaded  by  the  parabolic  load  were  suppressed.  This  suppres¬ 
sion  invokes  a  linear  edge  displacement  under  the  load. 
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a)  Two  Elements 


Loaded  Edge  Midpoint  Nodes 
are  Suppressed  Only  for 
Parabolic  Loading  Case 


b)  Eight  Elements 


nmm 

tmnm 


c )  Thirty-two  Elements 


Figure  113.  Idealization 


The  results  obtained  from  this  set  of  convergence  studies  are  presented 
in  Figures  114  and  115.  Figure  114  is  a  plot  of  the  membrane  displacement,  u^ 
at  the  middle  of  the  plate's  edge  versus  degrees -of -freedom  employed  in  the 
analyses.  The  reference  solution  (Reference  44)  is  designated  by  the  solid  line. 

Figure  115  presents  a  curve  of  the  membrane  displacement,  ux  and  stress  resul¬ 
tant,  Nx,  versus  the  edge  span  of  the  plate  for  the  idealization  shown  in  Figure  113  (c) 
(32  element  solution) .  It  should  be  noted  that  the  reference  solution  is  again  desig¬ 
nated  by  the  solid  lines  and  no  discernible  difference  between  solutions  can  be  detected. 

2.  Membrane  Gridwork  Influences 

The  second  illustration  which  utilizes  the  triangular  thin  shell  element  in 
a  structural  evaluation  will  be  the  following.  Again,  consider  a  thin  isotropic  square 
plate  loaded  with  a  self  equilibrating  parabolic  membrane  load  as  shown  in  Figure  112. 
This  illustration  will  involve  the  effect  that  the  shape  of  the  elements  used  in  the  struc¬ 
tural  idealization  has  on  the  determination  of  center  edge  displacement  of  the  mem¬ 
brane. 

The  six  idealizations  used  for  the  shape  study  are  shown  in  Figure  116.  Due  to 
conditions  of  symmetry  only  one  qi  .adrant  of  the  plate  was  analyzed.  The  midside 
nodes  which  were  loaded  by  the  parabolic  membrane  load  were  suppressed  in  this 
solution. 

The  results  obtained  from  the  subject  shape  studies  are  shown  in  Figure  117. 
These  solutions  indicate  that  the  displacement  values,  u  ,  obtained  for  the  middle 
of  the  plate's  edge  are  fairly  insensitive  to  the  shape  of  the  element  for  this  class  of 
problem. 


3.  Plate  Stress  Analysis 

The  third  illustration  which  utilizes  the  triangular  thin  shell  element  in 
a  structural  evaluation  will  be  the  following.  A  simply  supported  square  isotropic 
plate  with  a  uniform  normal  pressure  load  of  one  psi  is  shown  in  Figure  118  along 
with  its  material  properties  and  pertinent  dimensions. 

The  idealization::  used  for  the  finite  element  analyses  are  shown  in  Figure 
113.  Note  that  no  node  points  are  suppressed  in  this  analysis.  Three  different  grid 
sizes  were  employed  in  this  evaluation.  Two  element,  8  element,  and  32  element 
solutions  were  obtained  in  order  to  evaluate  convergence  characteristics.  Due  to 
conditions  of  symmetry  only  one  quadrant  of  the  plate  was  analyzed. 

Figure  119  is  a  plot  of  the  center  transverse  displacement  at  the  center 
of  the  plate  versus  the  degrees -of-freedom  employed  in  the  analyses.  The  reference 
solution  (Reference 45)  is  designated  by  the  solid  line.  Figure  120  presents  a  curve 
of  the  transverse  displacement,  wx  and  bending  moment,  Mx,  versus  the  center  span 
of  the  plate  for  the  idealization  shown  in  Figure  113  (c)  (32  element  solution) .  Again, 
it  should  be  noted  that  the  reference  solution  is  designated  by  the  solid  lines  and  no 
discernible  differences  between  solutions  can  be  detected. 
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Stress  Resultant  Nx  x  10“  (#/ in.) 


•  •  32  Element  Soln 

Exact  Soln 


Figure  115  Membrane  Displacement  And  Stress  Behavior 
Versus  Plate  Edge  Span 
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IDEALIZATION 

Figure  117  Shape  Study  Of  Triangular  Element  Membrane,  Square,  Isotropic, 
Self  Equilibrating  Parabolic  Edge  Loads 


®  Magic  Triangle 
■■■■  1  —  Reference  Alternate  Solution 

Figure  119  Convergence  Of  Triangular-Element  Plate,  Square 
Isotropic,  Simple  Supports,  Uniform  Pressure 


Legend : 

®  ^  Magic  Triangle 

Reference  Alternate  Solution 

Figure  120  Behavior /Triangular  Element  Plate,  Square,  Isotropic 
Simple  Support  Unit  Uniform  Load 
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4.  Plate  Gridwork  Influences 


The  fourth  illustration  which  utilizes  the  triangular  thin  shell  element  in 
a  structural  evaluation  will  be  the  following.  A  simply  supported  square  isotropic 
plate  with  a  uniform  normal  pressure  load  of  one  psi  is  shown  in  Figure  118  along 
with  its  material  properties  and  pertinent  dimensions. 

This  illustration  will  involve  the  effect  that  the  shape  of  the  elements  used 
in  the  structural  idealization  has  on  the  determination  of  the  maximum  displacement 
of  the  plate. 

The  six  idealizations  used  for  the  shape-  study  are  shown  in  Figure.116.  Due 
to  conditions  of  symmetry  only  one  quadrant  of  the  plate  was  analyzed. 

The  results  obtained  from  the  subject  shape  studies  are  shown  in  Figure 
121.  These  solutions  indicate  that  the  determination  of  the  plate's  center  transverse 
displacement  is  fairly  insensitive  to  the  shape  of  the  element  for  this  class  of  problem. 
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10.  DISCUSSION  AND  CONCLUSIONS 


A.  DISCUSSION 

Integrated  general  purpose  analysis  capabilities  of  the  MAGIC  System  class  sig¬ 
nal  a  major  advance  in  the  state-of-the-art  of  automated  tools  for  analysis.  The 
superior  cost  effectiveness  of  such  systems  over  conventional  multiple  special  purpose 
program  capabilities  is  compelling. 

This  assertion  of  superior  performance  from  large  scale  program  systems  may 
well  contradict  conclusions  drawn  from  experience.  Complexity  and  inefficiency  have 
long  been  concomitant  with  large  size  and  versatility  in  computer  programs.  Indeed, 
the  elimination  of  these  depreciating  effects  was  prerequisite  to  realization  of  the 
favorable  cost  effectiveness  of  the  MAGIC  System. 

Large  size  and  versatility,  without  excessive  complexity,  are  assumed  intrinsic 
to  the  MAGIC  System  in  subsequent  paragraphs,  as  attention  is  focused  upon  the  rela¬ 
tive  efficiencies  of  integrated  general  purpose  analysis  capabilities  and  multiple 
special  purpose  computer  program  analysis  capabilities.  This  is  to  presume  the  pre¬ 
requisite  elimination  of  the  greater  hindrance;  namely,  the  excessive  complexity 
which  choked  off  many  early  general  purpose  program  developments.  This  problem¬ 
atical  complexity  was  encountered  when  programs  of  simple  organization  grew  to 
press  upon  the  limits  of  computer  software  and  hardware  capabilities.  Extensions 
beyond  this  point  were  accomplished  by  intricately  coordinated  multiple  usage  of 
valuable  names  and  locations,  special  program  versions  with  omitted  features  and 
other  actions  which  accumulated  to  entangle  the  logic  and  data  storage  until  further 
modification  became  impractical. 

In  the  face  of  this  situation  increasingly  powerful  analytical  models  and  solution 
methods  were  formulated  and  numerical  implementation  demanded.  And,  as  is  often 
the  case,  sufficient  pressure  was  built  up  to  bring  about  the  technological  advances 
needed  in  the  computer  technologies. 

Advances  were  forthcoming  in  programming  technology  which  established  the 
technical  feasibility  of  a  truly  general  purpose  computer  program  system.  Advances 
in  computer  hardware  insured  the  economic  feasibility  as  the  technical  feasibility  was 
established  through  a  number  of  contributing  developments.  The  collective  result  of 
these  latter  deve’  ments  is,  in  a  word,  "organization."  Among  those  organizational 
characteristics  'atures  considered  essential  are,  the  breakdown  into  single  func¬ 
tion  modulei  .  ■<*  ,  igram  library  concept,  the  matrix  interpretive  system,  the 
SUBSYS  roui  .  In-depth  discussion  was  given  to  these  considerations  in  Section 

2  and  is  not  repeated  here.  Rather,  attention  is  given  to  the  benefits  which  accrue 
from  their  fulfillment  in  the  MAGIC  System. 
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It  is  appropriate  to  emphasize  at  this  point,  that  the  MAGIC  System  for  structural 
analysis  is  more  than  a  discrete  element  computer  program.  It  is,  in  one  sense,  a 
Problem  Oriented  Language  (POL)  which  enables  various  Analyst  specified  computa¬ 
tional  procedures.  And,  at  the  same  time,  it  is  designed  with  attendant  structural 
analysis  practices  evolved  from  applications  experience.  These  practices  are  dis¬ 
cussed  in  detail  in  subsequent  paragraphs.  The  point  of  interest  here  is  that  the 
efficiency  of  the  MAGIC  System  is  an  overall  efficiency  governed  more  by  men  than 
machines. 

The  more  comprehensive  the  comparison,  the  greater  the  advantage  shown  by 
the  integrated  general  purpose  analysis  capabilities  over  multiple  special  purpose 
program  capabilities.  In  nearly  all  cases  an  equitable  comparison  must  include  con¬ 
sideration  of  program  development  efforts  since  relevant  technologies  are  continuously 
advanced.  On  this  basis  the  integrated  approach  enjoys  the  greatest  relative  advantage. 
The  integrated  approach  is  also  superior  to  the  multiple  program  approach  when  con¬ 
sidering  only  factors  involved  in  utilization  of  operational  capability.  On  the  other 
hand,  shorter  execution  times  are  conceded  to  special  purpose  programs  without  dis¬ 
pute,  since  execution  efficiency  is  not  essential  to  the  case  for  the  greater  overall 
efficiency  of  integrated  analysis  capabilities. 

Attention  is  focused  now  on  the  impact  of  the  integrated  general  purpose  com¬ 
puter  program  approach  on  the  efficiency  of  the  many  processes  involved  in  mainte¬ 
nance  and  application  of  responsive  analysis  tools  in  support  of  a  broad  structural 
design  activity.  Program  maintenance  efforts  benefit  from  the  highly  modularized 
organizational  structure  to  an  even  greater  extent  than  the  initial  development  effort. 

In  the  initial  development,  functional  modules  are  established  against  the  require¬ 
ments  of  the  alternative  analysis  procedures  taken  collectively.  And,  since  an  exten¬ 
sive  commonality  exists,  multiple  repetitious  coding  is  avoided.  This  same  payoff  is 
derived  again  as  existing  modules  are  retired  in  favor  of  new  modules  which  offer 
improved  performance.  The  introduction  of  a  single  improved  module  is  reflected  to 
advantage  throughout  all  pertinent  analysis  procedures  of  the  computer  program  sys¬ 
tem.  The  option  exists  to  retain  alternative  modules  for  the  same  function  without 
sacrifice.  This  provides  useful  operational  flexibility  and  a  convenient  testbed  for 
various  candidate  procedures.  Alternative  procedures  can  be  evaluated  within  the 
system  without  disrupting  its  operational  status. 

The  foregoing  has  dealt  with  maintenance  of  existing  analysis  capability.  Main¬ 
tenance  is  also  interpretable  as  generalization  of,  and  addition  to,  the  overall  analysis 
capability.  Completely  new  analyses  can  be  implemented  with  the  addition  of  only 
those  functional  modules  absent  in  the  existing  capability.  For  example,  finite  element 
heat  conduction  analyses  are  possible  with  relatively  minor  modifications  to  the 
MAGIC  System. 
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The  benefits  derived  from  the  organization  of  a  general  purpose  computer  pro¬ 
gram  system  in  development,  maintenance,  generalization  and  extension  are  simul¬ 
taneously  important  disadvantages  associated  with  multiple  computer  program  analysis 
capabilities.  The  extensive  commonality  among  analyses  leads  in  this  latter  case  to 
the  repeated  development  of  coding  to  perform  a  given  function.  The  preparation  of 
special  versions  of  new  modules  and  the  introduction  of  these  into  a  multiplicity  of 
computer  programs  is  often  not  justified  and  the  overall  capability  is  depreciated. 

Another  particularly  important  handicap  borne  by  the  separate  programs  of  a 
multiprogram  capability  is  that  these  programs  can  not  command,  individually,  the 
provision  of  many  useful  special  features.  For  example,  useful  options  and  diagnostics 
are  usually  omitted  from  these  special  purpose  program  routines.  Also,  such  pro¬ 
grams  frequently  encounter  obstacles  such  as  machine  storage  capacity  which  must  be 
avoided  rather  than  surmounted  in  view  of  the  limited  applicability  of  the  program. 
Advancements  in  computer  software  and  hardware  are  further  considerations  of  im¬ 
portance  in  the  maintenance  of  an  analysis  capability.  Ti'ese  advancements  place 
multiple  program  capabilities  in  special  peril.  Those  programs  not  being  actively 
utilized  at  the  time  of  transition  in  software  or  hardware  are  easily  overlooked  and  in 
this  way  are  lost  from  the  overall  analysis  capability. 

No  single  factor  is  more  important  in  the  provision  of  a  responsive  analysis 
capability  than  documentation.  Engineering  documentation  must  delineate  analysis 
procedure,  input  data  and  output  data.  Programming  documentation  must  provide  for 
operation  and  modification  of  the  program. 

Consolidation  of  the  analysis  capability  into  a  general  purpose  program  results 
in  a  corresponding  favorable  consolidation  of  documentation.  Not  only  is  volume  re¬ 
duced  but  the  total  capability  is  described  uniformly  as  a  whole.  Small  programs  tend 
to  be  the  personal  tool  of  the  Initiator.  As  a  consequence,  the  documentation  prepared 
is  generally  inadequate  to  enable  general  usage.  This  situation  leads  to  extensive 
tutorial  instruction  to  realize  the  benefits  of  the  program  development.  At  the  very 
least,  multiple  program  capabilities  place  the  burden  of  assimilating  the  overall 
analysis  capability  from  the  individual  manuals  upon  the  user. 

The  foregoing  has  pointed  out  decisive  advantages  of  general  purpose  program 
systems  in  the  context  of  development  and  maintenance  of  analysis  capability.  The 
most  compelling  advantages,  however,  are  found  in  operation.  The  greater  efficiency 
of  the  MAGIC  System  relative  to  multiprogram  capabilities  for  analysis  stems  in 
large  measure  from  the  extent  of  the  analysis  process  which  is  covered.  Time  con¬ 
suming,  error  prone,  manual  transfers  of  data  between  special  purpose  or  single  step 
computer  programs  are  avoided.  The  integration  of  heat  conduction  and  thermal  stress 
analysis  within  a  single  system  can  circumvent  the  laborious  preparation  of  tempera¬ 
ture  data.  The  integration  of  stiffness  and  vibration  analyses  can  similarly  circumvent 
the  manual  transfer  of  stiffness  and  mass  data.  These  eliminations  of  manual  effort 
yield  reductions  in  calendar  time  which  Is  often  the  paramount  consideration  for  con¬ 
tribution  of  analysis  to  design.  This  is  not  to  say  that  long  continuous  executions  are 
desirable.  Execution  interruptions  enter  importantly  into  proper  utilization  of  the 
MAGIC  System. 
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The  MAGIC  System  is  designed  to  facilitate  good  structural  analysis  practices 
in  support  of  the  overall  structural  design  process.  Individual  design  organizations 
are  best  served  by  structural  analysis  practices  and  program  versions  which  are,  to 
some  degree,  distinct.  On  the  other  hand,  the  extensive  commonality  which  does  exist 
among  design  organizations  provides  strong  motivation  for  reviewing  the  effective 
structural  analysis  practices  and  supplemented  program  version  which  have  evolved  at 
Bell  Aerosystems. 

The  structural  analysis  process  begins  with  the  idealization  of  the  structure  into 
an  assemblage  of  finite  elements.  This  is  a  multistep  operation  if  the  structure  is 
first  separated  into  substructures.  Generally,  the  separation  into  substructures  is 
governed  by  the  physical  interconnections  of  the  major  structural  components.  The 
idealization  into  finite  elements  is  governed  by  variations  in  geometry,  dimensions, 
material,  applied  loading  and  boundary  conditions. 

Preprinted  input  data  forms  are  employed  to  simplify  and  thereby  improve  the 
reliability  of  the  input  data  specification.  These  preprinted  input  forms  associated 
with  the  MAGIC  System  are  an  important  improvement  over  card  image  forms  for 
frequent  as  well  as  infrequent  users  since  they  incorporate  automatic  data  generation 
features.  These  built-in  data  generation  features  are  supplemented  at  Bell  by  auxil¬ 
iary  (not  integrated  into  the  MAGIC  System)  data  generation  programs.  Some  of  these 
are  employed  routinely.  Others  are  extremely  simple  programs  written  for  a  single, 
problem  related  calculation.  Such  auxiliary  programs  are  frequently  employed  to 
advantage  in  the  generation  of  gridpoint  coordinates  with  reference  to  the  global  rec¬ 
tangular  coordinate  axes,  since  expression  of  these  can  require  extensive  tedious  cal¬ 
culation.  This  gridpoint  coordinate  data  set  should  be  interpreted  here  to  include 
points  for  specification  of  gridpoint  axes  transformations  and  stress  and  material 
angles  as  well  as  points  associated  with  degrees-of-freedom. 

The  first  MAGIC  System  execution  undertaken  is  to  confirm  the  assembled  input 
data  deck.  This  deck  is  read  and  the  implied  data  is  given  explicit  definition.  For 
example,  material  properties  are  extracted  from  the  Material  Library  and  gridpoint 
axes  transformations  are  generated  from  the  coordinate  table.  The  completed  data 
set  is  examined  in  this  preprocessing  execution.  All  data  items  are  stored  for  execu¬ 
tion  restart  and  printed  for  further  checking  by  the  analyst. 

The  preprocessing  execution  is  supplemented  at  Bell  to  include  the  generation  of 
a  magnetic  tape  which,  in  turn,  generates  a  plot  of  the  structural  model  on  an  automatic 
plotting  machine.  This  plot  enables  efficient  and  reliable  confirmation  of  the  two  most 
problematical  data  items;  namely,  the  gridpoint  positions  and  the  finite  element  con¬ 
nection  arrangement.  Beyond  this  point  the  structure  plot  is  a  useful  identifying  title 
sheet  for  the  printed  problem  output. 

The  next  phase  of  the  analysis  process  proceeds  via  a  restart  through  the  gener¬ 
ation  of  the  structural  matrices  for  stiffness,  stress,  loads,  assembly,  boundary  con¬ 
ditions,  etc.  Built-in  features  control  this  matrix  generation  to  selectively  form  only 
those  matrices  required  for  the  current  analysis.  Completion  of  the  matrix  generation 
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phase  signals  exit  from  the  Structural  System  Monitor.  This  is  an  interface  point  be¬ 
tween  matrix  abstraction  instruction  statements,  and,  therefore,  a  point  for  optional 
interruption  of  the  execution  to  examine  the  system  level  matrices.  This  interruption 
is  used  only  infrequently  at  Bell. 

Calculation  proceeds  under  the  FOilMAT  System  to  the  governing  matrix  equa¬ 
tion  and  thence  to  the  solution  for  the  displacement  vectors  for  all  load  conditions. 

For  some  problems  execution  may  be  terminated  at  this  point.  For  many  other 
problems  the  validity  of  the  analysis  can  be  assessed  against  these  displacement 
results  and  an  execution  interruption  is  justified  by  the  computational  invest¬ 
ment  required  for  the  secondary  results.  Ideally,  the  deformed  structure  should 
be  plotted  to  facilitate  interpretation  of  the  predicted  displacement  behavior. 

The  analysis  proceeds  from  the  displacement  solution,  with  or  without  interrup¬ 
tion,  to  calculation  and  print  of  the  remainder  of  the  output  data  items;  namely,  reac¬ 
tions,  forces,  stresses,  etc.  This  is  the  conventional  point  of  termination  of  finite 
element  analyses.  However,  a  number  of  relatively  simple  auxiliary  programs  are 
used  to  advantage  at  Bell  to  relieve  the  burden  this  output  places  on  the  stress  engi¬ 
neers.  As  in  the  case  of  the  input  data  generation  auxiliary  programs,  some  of  the  auxi¬ 
liary  output  data  reduction  programs  are  employed  repeatedly  and  others  are  special 
to  a  single  problem.  The  functions  of  these  programs  include  such  things  as  principal 
stress  calculations  and  margin  of  safety  determinations.  Auxiliary  programs  which  do 
nothing  but  selectively  print  and  label  output  data  items  are  also  helpful  for  large 
problems. 

Several  comments  on  the  evaluation  of  output  data  are  warranted  in  concluding 
discussion  of  good  structural  analysis  practices.  The  examination  of  output  by  the 
Analyst  should  be  initiated  under  the  presumption  that  an  error  exists  with  confidence 
in  the  validity  of  the  analysis  accumulating  as  the  examination  proceeds.  Given  a  com¬ 
plete  set  of  output,  attention  should  first  be  given  to  the  gridpoint  force  balances  and 
reactions.  Assured  that  no  unintended  reactions  exist  and  that  residuals  are  negligibly 
small,  the  displacement  states  should  be  examined.  If  the  general  deformed  configura¬ 
tion  does  not  expose  any  inconsistencies,  confirmation  is  completed  by  examination  of 
the  more  extensive  presentation  of  force  and  stress  data. 

The  foregoing  discussion  has  focused  upon  development,  maintenance  and  utiliza¬ 
tion  considerations  important  to  the  favorable  cost  effectiveness  of  the  present  MAGIC 
System  for  structural  analysis.  Further  evolution  of  this  system  can  be  expected 
which  will  continue  to  improve  its  relative  advantage.  Updated  versions  of  the  MAGIC 
System  will  be  compatible  with  all  features  developed  in  connection  with  prior  versions. 

B.  CONCLUSIONS 

It  is  concluded  that  the  subject  MAGIC  System  provides  a  capability  for  structural 
analysis  equivalent  to  that  of  the  predecessor  programs  delivered  under  Contract 
AF33(657)8963,  taken  collectively.  The  satisfactory  achievement  of  this  overall  objec¬ 
tive  is  given  substantiation  by  a  number  of  subsidiary  conclusions.  Specifically,  it  is 
concluded  that: 
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(1)  The  finite  element  library  enables  effective  idealization  of  most  linear 
structures. 

(2)  Computational  procedures  attendant  to  the  MAGIC  System  enable  the  conduct 
of  linear  displacement  and  stress  analyses  in  the  presence  of  general  pre¬ 
strain  and  thermal  loading  as  well  as  distributed  and  concentrated  mechani¬ 
cal  loading. 

(3)  The  stability  analysis  procedure  provided  in  the  MAGIC  System  enables  the 
prediction  of  critical  load  levels  for  general  framed  structures. 

(4)  The  preprinted  input  data  forms  facilitate  the  rapid  and  reliable  specifica¬ 
tion  of  problem  data. 

(5)  The  computer  program  organization  of  the  MAGIC  System  effectively 
utilizes  the  FORMAT  System  and  is  well  suited  to  generalization. 
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